AD-A036  972  ARMY  ENGINEER  WATERWAY*  EXPERIMENT  STATION  VICKSBURG  MISS  P/9  8/6 

LAKE  ERIE  INTERNATIONAL  JETPORT  MODEL  FEASIBILITY  INVESTIGATION— ETC <U) 
APR  77  D C RANEY » D L DURHAM#  H L BUTLER 
UNCLASSIFIED  WES-TR-H-7G-6  NL 


Prtptrtd  for  Lake  Erie  Regional  Transportation  Authority 
Cleveland,  Ohio  44113 

Undor  Task  17  of  LERTA  Second-Phase 
Airport  Feasibility  Study 


TECHNICAL  REPORT  H-74-6 


LAKE  ERIE  INTERNATIONAL  JETPORT 
MODEL  FEASIBILITY  INVESTIGATION 

Report  17-4 

NUMERICAL  MODEL  FEASIBILITY  STUDY 


Donald  C.  Raney,  Donald  L.  Durham,  H.  Lee  Butler 
Hydraulics  Laboratory 

U.  S.  Army  Engineer  Waterways  Experiment  Station 
P.  O.  Box  631,  Vicksburg,  Miss.  39180 


April  1977 
Report  4 of  a Series 


‘‘A 


Approved  For  Public  Release;  Distribution  Unlimited 


Destroy  this  report  when  no  longer  needed.  Do  not  return 
it  to  the  originator. 


1/ 

The  preparation  of  this  document  was  financed  in  part  through  an  Airport  Master 
Planning  Grant  from  the  Department  of  Transportation,  Federal  Aviation  Administration, 
under  the  Planning  Grant  Program  as  provided  in  the  Airport  and  Airway  Development 
Act  of  1970.  The  contents  of  this  report  reflect  the  views  of  the  6p»e  Western  Hi  inn 
•ye<HPrty,  which  is  responsible  for  the  facts  and  accuracy  of  the  data  presented  herein, 
i The  contents  do  not  necessarily  reflect  the  official  views  or  policy  of  the  FAA, 

« tf 

Department  of  the  Army,  Corps  of  Engineers,  or  Environmental  Protection  Agency. 
Acceptance  of  this  report  by  the  FAA  does  not  in  any  way  constitute  a commitment 
on  the  part  of  the  United  States  to  participate  in  any  development  depicted  therein, 
nor  doe*,  it  indicate  that  the  proposed  development  is  environmentally  acceptable  in 
i'  accordance  with  Public  Laws  91-190,  91-258,  and/or  90-495. 


DEPARTMENT  OF  TME  ARMY 

WATERWAYS  EXPERIMENT  STATION.  CORPS  OF  ENGINEERS 
P O.  BOX  631 

VICKSBURG.  MISSISSIPPI  3S1BO 

.v  Hire*  TO.  WESAR  27  April  1977 

Errata  Sheet 
No.  1 

LAKE  ERIE  INTERNATIONAL  JETPORT  MODEL  FEASIBILITY  INVESTIGATION 
NUMERICAL  MODEL  FEASIBILITY  STUDY 

Technical  Report  H-74-6 
Report  17-4 

April  1977 

On  the  inside  of  the  cover,  in  line  4 of  the  disclaimer,  "Case  Western 
Reserve  University"  should  read  "U.  S.  Army  Engineer  Waterways  Experi- 
ment Station." 


SECURITY  CL  ASSI  F|C  ATlON  OF  THIS  PAGE  f'WTiwn  Pete  f-'ntered) 


L 


l7-_AHIHQftUJ 


REPORT  DOCUMENTATION  PAGE 


REPORT  number 

: rt  " -U  ^ 


2.  GOVT  ACCESSION  NO 


4 T|  T L F < and  Subtitle) 


LAKE  l K Jim  RNATIONALJETPf  *PT  ^ODEL  j’EASIBI I, TTY  I 

- .-r:  ;at5  numer]  :ai  m dei  feasibility  study  I 

- t?T  <*  T_  'V  / 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


RECIPIENT'S  : ATALOG  NUMBER 


S TYPE  OF  REPORT  A PERIOD  COVERED 

Report  )i  of  a series  Z-'" 


i PERFQRMINOrORO-  REPORT  NUMBER 

j 'rechni 


ical  /Repyt  :Ji— 7 U— G 

B CONTRACTOR  NUMBERCa) 


Donald  C.  Raney. 

Donald  L.  Durhhrc 
H.  Lee  Butler 

V PFftt'oRMING  ORGANI  Z ATlON  NAME  ANd’aDDRESS 

U.  S.  Army  Engineer  Waterways  Experiment  Station 

Hydraulics  Laboratory 

P.  L.  Box  631,  Vicksburg,  Miss.  391 -'>0 


WE  - - L iilAzrtdl 

’5S  \ '0.  PROGRAM  BREMEN  T.  PROJECT,  TASK 


11.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Lake  Erie  Regional  Transportation  Authority 
33  Public  Square,  Suite  1015 
Cleveland,  Ohio  UU.ll 3 


T5  MONITORING  AGENCY  NAME  ft  ADDRESS^//  different  from  Controlling  Office) 


AREA  A WORK  UNIT  NUMBERS 
/:  /• 


12.  REPORT  DATE 

n 1 vnfrrmaaj  T7 

( ' J'n.-i tg 


HUMBER11?) V pages 

151 


15.  SECURITY  CLASS,  (of  thle  report) 

Unclass i f ied 


15*  DECLASSIFICATION  DOWNGRADING 
SCHEDULE 


16  DIS-RI0U  TION  STATEMENT  (of  t hla  Report) 


Approved  for  public  release;  distribution  unlimited. 


17.  DISTRIBUTION  STATEMENT  (of  the  mbetract  entered  In  Block  30,  If  different  from  Report) 


IS  SUPPLEMENTARY  notes 


19  KEY  WORDS  ( Continue  on  revetee  aide  If  neceaemry  end  Identify  by  block  number) 

Airports 
Hy dr  dynamics 
Lake  Erie 

Mathematical  models 


20  ABSTRACT  ("Continue  a n rover ee  aide  It  nocoeeory  and  Identify  by  block  number) 

An  integral  part  of  the  feasibility  assessment  of  a proposed  offshore 
Jet]  rl  Lti  need  Cleveland,  hio,  is  the  investigation  >f  the  hydrodynamics  f 
Lake  Erie  to  aid  in  determining  the  effects  of  the  structure  on  sucli  phenomena  as 
ichii  ■ , 1 ra  , and  Lake  circulation.  as  1st  Ln  leti  m Lning  these 

eff<  :ti  , thi  feasibility  of  using  numer  teal  modeling  techniques  tra  invest i ;at<  i . 
Numerical  models  that  appeared  capable  of  predicting  the  extent  and  magnitude  of 
hydr  iynamic  chai  -•  produced  by  the  proposed  Jet  port  wen  reviewed.  Basi  1 u]  >n« 


Af  X 


1 


DO  ,'r*n  1473 


EDITION  or  I MOV  6»  IS  OBSOLF  TE 


Unclasp i fied 


SECURITY  CL  ASSl  F|C  AT]ON  OF  THIS  PAGf  When  Pete  Entered ) 


■ / 


l no  I ; V i » • i 

SECURITY  CLASSIFICATION  of  This  PAGEfW7»«n  Dml»  Rnffd) 

20.  ABSTRACT  (Continued). 

this  investigation,  a determination  was  made  f the  fea  . t Llity  of  apj 
numerical  models  to  the  pi oblems  of  sei chine,  storm  surge,  and  lake  circulation 
in  Lake  Erie.  The  numerical  model  feasibility  study  was  restri :ted  1 
consideration  of  existing  state-of-the-art  models.  No  extensive  model  lev<  . ; - 
ment  was  undertaken  in  this  feasibility  study;  however,  areas  where  additional 
numerical  development  is  required  were  identified.  The  theoretical  Limitat 1 • , 
verification,  resolution  capabilities,  and  accuracy  of  the  model  as  w^n  as  th< 
cost  of  application  of  the  model  wen-  considered. 


Unclass i fied 


SECURITY  CLASSIFICATION  OF  THIS  PAGEfWTiFn  D«r«  Fnr»«f) 


THE  CONTENTS  OF  THIS  REPORT  ARE  NOT  TO  BE 
USED  FOR  ADVERTISING,  PUBLICATION,  OR 
PROMOTIONAL  PURPOSES.  CITATION  OF  TRADE 
NAMES  DOES  NOT  CONSTITUTE  AN  OFFICIAL  EN- 
DORSEMENT OR  APPROVAL  OF  THE  USE  OF  SUCH 
COMMERCIAL  PRODUCTS. 


I 


PREFACE 


This  study  was  sponsored  by  the  Lake  Erie  Regional  Transportation 
Authority  (LERTA),  Cleveland,  Ohio,  as  a part  of  the  model  feasibility 
investigation  being  conducted  at  the  U.  S.  Army  Engineer  Waterways 
Experiment  Station  (WES).  The  WES  investigation.  Task  IT  of  the  LERTA 
investigation,  is  a portion  of  the  second-phase  airport  feasibility 
study  undertaken  by  LERTA  to  evaluate  proposed  airport  sites,  one  of 
which  is  in  Lake  Erie  near  Cleveland.  The  numerical  model  feasibility 
study  is  associated  with  the  selection  and  preliminary  design  of  the 
necessary  numerical  models  for  studying  various  phenomena  considered 
pertinent  to  an  offshore  jetport  site. 

This  report  was  prepared  by  Dr.  Donald  C.  Raney,  Dr.  Donald  L. 
Durham,  and  Mr.  H.  Lee  Butler  of  the  Wave  Dynamics  Division  (WDD), 
Hydraulics  Laboratory  (HL),  WES,  under  the  general  supervision  of 
Dr.  R.  W.  Whalin,  Chief,  WDD,  and  Mr.  H.  B.  Simmons,  Chief,  HL. 

Dr.  Raney  is  a Professor  in  the  Aerospace  and  Mechanical  Engineering 
Department,  University  of  Alabama,  and  was  assigned  to  WES  under  terms 
of  the  Intergovernmental  Personnel  Exchange  Act  during  conduct  of  this 
study  and  preparation  of  the  report. 

The  Directors  of  WES  during  the  conduct  of  this  investigation  and 
the  preparation  and  publication  of  this  report  were  COL  G.  H.  Hilt,  CE, 
and  COL  J.  L.  Cannon,  CE.  Technical  Director  was  Mr.  F.  R.  Brown. 
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CONVERSION  FACTORS,  METRIC  (Si)  TO  U.  S.  CUSTOMARY  AND 
U.  S.  CUSTOMARY  TO  METRIC  (Si)  UNITS  OF  MEASUREMENT 


Units  of  measurement  used  in  this  report  can  be  converted  as  follows: 


Multiply 

To  Obtain 

Metric 

(SI)  to  U.  S. 

Customary 

1 

metres 

3.280839 

feet 

kilometres 

0.6213711 

miles  (U.  S.  statute) 

metres  per  second 

3.280839 

feet  per  second 

square  centimetres 
per  second 

0.1550 

square  inches  per  second 

i'elvins  or  Celsius 

degrees 

9/5 

Fahrenheit  degrees* 

U.  S. 

Customary  to  Metric  (Si) 

inches 

2. 51* 

centimetres 

feet 

0.301*8 

metres 

miles  (U.  S.  statute) 

1.60931*!* 

kilometres 

square  feet 

0.0929030U 

square  metres 

cubic  feet 

0.02831685 

cubic  metres 

pounds  (force)  per 
inch  absolute 

square 

689!*.  757 

pascals 

feet  per  second 

0.30!*8 

metres  per  second 

miles  per  hour 
(U.  S.  statute) 

1.6093!*!* 

kilometres  per  hour 

degrees  (angle) 

0.0171*5329 

radians 

degrees  Fahrenheit 

5/9 

Celsius  degrees  or  kelvins** 

* To  obtain  Fahrenheit  (F)  temperature  readings  from  Celsius  (C)  read- 
ings, use  the  following  formula:  F = 9/5(C)  + 32.  To  obtain  Fahren- 
heit readings  from  Kelvins  (K),  use:  F = 9/5(K  - 273.15)  + 32. 

**  To  obtain  Celsius  (C)  temperature  readings  from  Fahrenheit  (F)  read- 
ings, use  the  following  formula:  C = (5/9)(F  - 32).  To  obtain  Kel- 
vin (K)  readings,  use:  K = ( 5/9 ) ( F - 32)  + 273.15. 
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J.AKK  ERIK  ] KTERNATIONAI.  .JETPORT 
MODEL  FEAUIMLJTY  INVEDTIOATION 


NUMERICAL  MODEL  FEASIBILITY  AND  SEICHE  STUDY 


PART  I:  INTRODUCTION 


Background 

1.  Tlie  Lake  Erie  Regional  Transportation  Authority  (LERTA)  is 
conducting  a feasibility  and  site  selection  study  for  a major  hub  air- 
port in  the  Cleveland  Service  Area.  One  of  the  sites  being  evaluated  is 
offshore  in  Lake  Erie  near  Cleveland,  Ohio.  As  a part  of  the  feasi- 
bility analysis  of  an  offshore  site,  the  U.  S.  Army  Engineer  Waterways 
Experiment  Station  (WES)  is  conducting  a model  feasibility  investigation. 
Objectives  of  the  WES  study'1'  are: 

a.  Compilation  of  available  data  on  wave  activity  (wind  waves, 
seiches,  and  tides)  and  mass  circulation  in  Lake  Erie  with 
particular  emphasis  on  effects  of  these  phenomena  in  and 
around  Cleveland. 

b.  Selection  and  preliminary  design  of  necessary  hydraulic 
models  for  studying  various  phenomena  considered  pertinent 
to  the  jetport  site. 

c_.  Evaluation  and  preliminary  application  of  analytical  and/or 
numerical  models  of  seiching  and  mass  circulation  in  a lake 
to  the  jetport  study. 

The  five  study  tasks  in  the  WES  investigation  are: 

a.  Synthesis  of  available  data  primarily  concerning  wave 
climate,  mass  circulation,  general  shoreline  character- 
istics, and  general  features  of  erosional  problems. 

b.  Lake  seiche  analysis. 

c_.  Wave  diffraction  and  refraction  analyses  including  quali- 
tative effects  on  shore  erosion. 

d_.  Analytical  mass  circulation  analysis. 

e_.  Preliminary  design  of  necessary  hydraulic  models. 

This  report  is  the  fourth  report  in  the  series  published  under  the 
general  title  "Lake  Erie  International  Jetport  Model  Feasibility 
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nv<  tigati  n. " result  f studj  ta  and  a port i f ;tudj 

i_  art  presented.  Inder  studj  task  1 >nlj  thi  initial  select! 

: ui  rica  models  t b«  applied  is  presented  in  this  rej  >rt . De- 
■ 1 results  fr  m applicatd  i f thesi  n am  rica]  n iels  will  b<  pre- 
sented in  later  report.'. 


lruHiTn  i '<  • r i n i ti-.n 

. . An  integral  part  of  the  feasibility  assessmi  nt  fa  pro]  . ■ : 
hore  jetport  site  near  Cleveland  Harbor  is  the  invest  igatioi 
hydrodynamics  of  Lake  Erie  to  aid  in  determining  effects  of  the  struc- 
ture on  such  phenomena  as  seiching,  storm  surge,  -nil  Lake  :ircul<  ti  in. 

Ti  assist  in  determining  these  effects,  the  feasibility  if  using  num<  r- 
;a]  ■■  :•  Ling  techni  jues  was  investigated.  Numer ica]  n iels  that 
appeared  capable  of'  predicting  the  • xtent  and  magnitude  of  hydr  lyn  i 
; iuced  by  thi  proposed  jetport  were  reviewed.  Based  i]  n 
this  investigation  of  existing  state-of-the-art  numeri  :a]  n iels,  at 
assessment  of  the  feasibility  of  applying  numerical  models  to  the 
pr  ; Lems  of  seiching,  storm  surge,  and  wind-driven  circulati  n it 
: r i ■ is  made.  Any  models  selected  for  this  purpose  must  pr  tuffi- 

• ' ■ nl  resolution  to  allow  the  hydrodynamic  changes  induced  bj  th<  jet- 
port  to  be  observable.  This  survey  and  evaluation  could  not  cover  all 
ting  it  lei  r models  under  development ; therefore,  vari  is  n I 
available  from  the  literature  and  personal  knowledge  of  tlv  authors 
were  considered  at  the  time  of  this  survey. 

3.  One  phenomenon  investigated  during  this  study  Ls  seiching. 
Seiches  are  long-period  oscillations  of  the  lake  surface  about  the 
• ■ evel.  Thesi  tending  waves  ar<  fo rmi id  after  a wind  bl ow ing 

across  the  Laki  subsides  and  t ; e setu]  of  the  water  surface  Is  n . ■ ■ r 
maintained  bj  the  wind  stress . Thi  free  iscillations  have  period 
arc  dependent  upon  the  horizontal  and  vertical  dimensions  of  the  Ink--, 
fricti  n,  and  the  number  of  nodes  of  the  stai  ling  wave,  thal  Ls,  Lines 
deviation  of  the  free  surface  from  Its  undisturbed  value  is  s«  r . 
Seiching  affects  both  water  level  and  mass  :irculat  i in  patti  ms.  Li  was 
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ressa  rj  • invest igati  th<  peri  is , mode  configurations , and  vel oc  ! ty 
limes  f r th<  fre<  ■ i ations  .ak<  ri(  l tl  w it!  at  3 withoul  th< 
jet i rl  structure.  his  investigate  i provides  relative  information  for 

• ; • : f th<  ■ ffects  of  such  a structure  on  the  mode  confi  rural  i ns 

and  peri  i f th<  fre<  ; illations  f Lake  Er Le. 

An  ther  phen  men  n s metim<  i ass  :iat<  3 with  seiching  is  si  nr 

surge.  Surge  is  i Lake  Level  fluctuati  t :ause  . bj  the  wi nd  stre 

acc  mpanying  a n ving  si  ru  ;ystem.  A wind  hlowin  ; v he  Lake  ex<  rl  ; 

a horizontal  force  on  the  lake  surface  and  induces  a surface  current  in 

the  general  direction  of  the  wind.  There  currents.  are  i t:i[  eded  in 

w-wat<  r ar  as,  thus  causing  the  lake  level  to  rise  down  wind  and 
the  windward  si  le.  Terms  used  to  describe  these  fluctuat 3 t 
el  up"  and  "set  lown."  The  Lake  surfac  fluctuati  ns  and  a :i- 

at'.-d  • irrents  produced  by  moving  storm  systems  are  considered  both  with 

i without  the  .jetport  structure. 

! . A third  phenomenon  of  importance  in  this  investigation  is  mass 
circulation  Ln  the  lake.  The  main  driving  force  producing  mass  :ircu- 

• ii  Lake  Erie  the  winds.  The  wind-generated  currents  can  be 
ite  1 ini  tw  *a1  - ries:  quasi-steady  and  time  dependent.  In  a 

sens  , the  re  are  no  steady-state  currents  in  the  lake  sii  ••  the 

:tor  :ausing  the  • irri  nti  are  usually  varying  in  time.  However, 

there  are  periods  when  the  wind  is  varying  sufficiently  slowly  to  permit 

teady-  tate  analysis  to  be  approximately  valid.  A steady~stat< 

numerical  model  based  on  an  extension  of  the  shallow-lake  theory  of 
2 

Welander  should  be  valid  during  the  fall  and  winter  months  since  the 
lake  is  well  mixed  during  this  period  due  to  wind-driven  turbulence  and 
thermal  convection. 

6.  When  these  well  mixed  conditions  are  not  appli  :al  L<  , as  Ln  th< 

urn iths,  a stratified  model  is  required  to  investigate  the  mass 

circulation.  A stratified  model  provides  for  a cooler  bottom  layer 
(hypolimnion)  separated  from  the  top  layer  (epilimnion)  by  a thin 
thermocline.  Of  particular  interest  is  any  effect  of  the  jetport 
structure  on  location  of  the  thermocline  and  the  effect  on  hypolimnion 
and  epilimnion  circulation  produced  by  typical  summer  wind  conditions. 
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7.  The  numerical  model  feasibility  study  was  restricted  to  the 
consideration  of  existing  state-of-the-art  models.  No  extensive  model 
development  wa.;  undertaken  in  this  feasibility  study;  however,  areas 
where  additional  numerical  development  are  required  were  identified. 

The  theoretical  limitations,  veri f ication , and  accuracy  of  the  model  as 
well  as  the  cost  of  application  of  the  model  were  considered. 
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PART  II:  SEICHE  ANALYSIS 


8.  In  this  part  of  the  report,  an  analytical  and  numerical  study 
of  the  effects  of  a jetport  in  Lake  Erie  near  Cleveland,  Ohio,  on 
seiches  in  the  lake  is  presented.  For  the  first  five  modes  of  free 
oscillations  in  Lake  Erie,  the  mode  configuration,  period,  and  velocity 
regimes  both  with  and  without  a jetport  are  calculated.  These  data  are 
used  as  relative  information  for  estimating  the  effects  of  a jetport 
structure  on  seiches  in  Lake  Erie. 

Mathematical  Formulation 

9.  As  previously  defined,  seiches  are  vertical  oscillations  of  the 
water  surface  about  the  mean  lake  level.  For  a body  of  water  with 
variable  depth  and  with  a maximum  depth  which  is  small  compared  with  the 
wavelength  of  the  oscillation,  a boundary  value  problem  for  seiches  in 
Lake  Erie  can  be  formulated  from  the  long-wave  equation: 


where 

n(x,y,t)  = surface  displacement  as  a function  of  time  and  space* 
t = time 

g = acceleration  due  to  gravity 
b(x,y)  = water  depth 
x,y,z  = Cartesian  coordinates 

The  coordinate  system  is  shown  in  Figure  1 in  which  the  undisturbed  free 
surface  is  the  plane  z = 0 , where  z is  the  vertical  distance,  the 
bottom  is  z = -h(x,y)  , and  the  surface  displacement  is  z = n(x,y,t)  . 
For  free  oscillation  in  a basin,  the  surface  displacement  can  be  assumed 
to  be  harmonic  in  time. 

* For  convenience,  symbols  ard  unusual  abbreviations  are  listed  and 
defined  in  the  Notation  (Appendix  A). 
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Figure  1.  Coordinate  system  and  surface  profile 


Thus,  n can  be  represented  as: 

n(x,y,t)  = Re[s(x,y)  exp  (-iwtj]  (2) 

where 

Re  = real  number 
£ = wave  amplitude 
i = imaginary  number 

w = 2m/T  where  T = period  of  oscillation 
Based  on  Equations  1 and  2,  the  governing  partial  differential  equation 
for  free  harmonic  oscillations  of  a body  of  water  becomes 


3_ 

3x 


h(x,y) 


T 3£(x,y )" 

] + 1_  r 

L 3x 

h(x,y 


+ 


A2£(x,y) 


= 0 


(3) 


p p 

where  X is  an  eigenvalue  expressed  as  X = w /g.  The  boundary  condi- 
tion for  this  problem  is  that  the  normal  component  of  velocity  at  the 
basin  boundary  be  equal  to  zero.  Since  Equation  3 (special  form  of  the 
Helmholtz  equation)  is  the  equation  for  a standing  wave,  this  boundary 
condition  can  be  expressed  as 


K = 

3n 


0 


(M 
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where  n = unit  normal  to  boundary.  Therefore,  Equations  3 and  h are 
the  partial  differential  equation  and  the  boundary  condition,  respec- 
tively, for  tiie  boundary  value  problem  of  the  free  oscillations  of  a 
body  of  water. ^ The  x and  y components  of  velocity  (u  and  v , 
respectively)  can  be  obtained  by  integrating  over  a suitable  time 
interval  the  equations  of  x and  y momentum  with  the  pressure  gradient 
expressed  as  the  gradient  of  the  surface  displacement.  These  components 
can  be  shown  to  have  the  following  form: 


Thus,  a solution  to  the  boundary  value  problem  for  the  surface  displace- 
ment or  the  gradient  of  the  surface  displacement  having  been  obtained, 
the  velocity  components  can  be  calculated  from  Equation  5. 

Numerical  Scheme 

Finite  element  method 

10.  The  boundary  value  problem  of  Equations  3 and  h can  be  solved 
analytically  for  basins  of  simple  geometry  and  topography.  However,  for 
basins  of  complex  geometry  and  variable  topography  such  as  Lake  Erie, 
solutions  to  this  problem  must  be  obtained  numerically.  The  numerical 
approach  chosen  to  solve  this  boundary  value  problem  as  applied  to  Lake 
Erie  is  a variational  approach  using  the  finite  element  method  (FEM).^ 
The  FEM  is  a discrete  approximation  procedure  applicable  whenever  a 
variational  principle  can  be  formed. J This  principle  is  expressed  as 

an  integral  over  the  region  under  consideration  which  is  discretized 
into  small  regions  known  as  finite  elements.  Within  each  element  the 
dependent  variable  (surface  displacement)  is  approximated  by  a local 
Taylor  series  expansion.  Since  there  are  no  constraints  imposed  on  the 
shape  of  each  element,  the  FEM  can  conveniently  and  easily  solve  bound- 
ary value  problems  involving  irregularly  shaped  boundaries. 

11.  The  variational  formulation  for  the  boundary  value  problem  of 
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s 


the  free  oscillations  of  a body  of  water  consists  of  the  functional 
integral 


where  A is  the  surface  area  of  Lake  Erie.  For  a stationary  value  of 
X (maximum  or  minimum)  with  respect  to  £(x,y)  , the  Euler-Lagrange 
condition  of  the  functional  integral  is  identical  with  Equation  3-  The 
identical  boundary  condition  for  Equation  I4  is 


h 


(7) 


12.  The  solution  of  £(x,y)  for  Equation  6 is  found  by  applying 

the  finite  element  techniques  developed  for  what  are  known  as  constant 

strain  triangles  (CST).  The  FEM  coding  for  the  program  was  developed 

using  linear  element  geometries  and  interpolation  functions  that  de- 

L 

scribe  CST  as  outlined  by  Desai  and  Abel.  Dividing  A into  N sub- 

0 

regions  or  elements,  A where  e = 1 ...n  , the  shape  of  an  element 
can  be  defined  by  a number  of  nodes  (vertices)  that  connect  the  element 
to  other  elements  of  the  grid.  Equation  6 can  now  be  expressed  at  the 
element  level  as 


V -V 

ff  i HeP 

1 

+ he 

9£U,y) 

x ~ 2->~ 

IJ  2lh  L 

9x  J 

+ n 

. 9y 

e=l 

A 

where  he 

= the  mean 

depth  for 

each 

element . 

t 2 


- xV( 


x,y)j 


dx  dy  ( 8 ) 


element,  £(x,y)  can  be  assumed  to  be  approximated  by  a Taylor  series 
expansion : 


£(x,y)  = a]  + a^x  + a^y 


(9) 
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where  , a,  , arid  a ^ are  unknown  coefficients  of  the  expansion. 
By  defining  the  matrix  vector  (£}e  as 


and  the  coordinates  of  the  three  nodes  forming  the  elements  as  (x.,y.), 

(x  ,y  . ) , and  (x  ,y  ) where  i,  j,  and  k are  indices  of  these  nodes, 
j .1  k k 

Equation  9 can  then  he  expressed  in  matrix  form  for  each  element  as: 


lx.y . 

a-, 

i i 

1 

lxiyj 

a2 

_lxkyk 

_a3_ 

By  expressing  the  a matrix  in  terms  of  the  node  coordinates,  a matrix 
[N]  , which  is  known  as  the  "interpolation"  or  "shape"  function,  can  be 
defined  as 


[»i  - [vYVl 


where 


'i-[' 


Vj  - Vk1  * (yk  - 1) 

e 

of  the  element  A 


)x  + (x.  - x^)yJ^2A  and  A = area 


N.,N  = similar  expressions  obtained  by  the  cyclical  permutation 
J k 

of  i , j , and  k 

Hence  £(x,y)  (Equation  9)  can  be  expressed  for  any  point  within  the 
element  as 


S(x,y)  = [N](Oe 


(10) 


Now,  with  C uniquely  and  continously  defined  throughout  the  region  of 
0 

A , the  functional  x Cfin  be  minimized  or  maximized  with  respect  to 
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the  nodal  values  of  £ . Substitution  of  Equation  10  into  Equation  8 
and  setting  the  first  variation  of  x with  respect  to  any  nodal  point 
m equal  to  zero  gives 


for  m = 1...M  , where  M is  the  total  number  of  node  points.  Equa- 
tion 11  represents  a system  of  M linear  algebraic  equations  in  terms 
of  the  M values  of  £ of  the  nodal  points,  in  which  the  values  of 
Cm  form  the  solution  of  Equation  3.  The  system  of  equations  formed  by 
Equation  11  can  be  expressed  in  general  matrix  format  for  the  entire 
region  of  interest  as 


[K]  (C)  = A2[MK0 


(12) 


where  X is  the  eigenvalue  expressed  as 

-.2  _ 2, 

X - It!  /g 

and  [ K ] = element  stiffness  matrix 
[M]  = lumped  mass  matrix 

Details  of  the  expression  of  individual  terms  of  [K]  and  [M]  matrices 
can  be  found  in  Reference  7.  After  assemblage  of  Equation  12,  the 
boundary  condition  (Equation  7)  is  inserted  at  the  appropriate  nodal 
values . 

13.  Equation  12  is  a form  of  an  eigenvalue  problem.  Thus,  for 

real,  square  matrices  of  order  M there  will  be  n independent 

real  eigenvalues,  X^  , where  n = 1...M  , which  satisfy  Equation  12 

and  define  the  natural  frequencies  of  the  system.  For  each  X there 

n 

will  be  an  eigenvector  or  "mode"  {£}  in  which  the  relative  magnitudes 
of  nodal  displacements  are  fixed  hut  not  their  absolute  values 

lh 


(solution  not  unique).  To  solve  this  eigenvalue  problem  by  direct  matrix 
procedures  requires  the  conversion  of  Equation  12  into  the  standard  form 

1 H ] { Z . } = X.{Z.)  , i = 1,2, ...M  (13) 


or 


[ H ] [ Z ] = [ A ] [ Z ] 

where 

[H]  = M x M square  matrix  of  coefficients 

[Z]  = modal  matrix,  where  columns  are  eigenvectors  { Z ^ } of  (II) 
[A]  = diagonal  matrix  of  eigenvalues 


Procedures  for  reducing  Equation  12  to  standard  form  (Equation  13)  are 
presented  in  detail  in  References  7 and  8. 

lk.  Once  {£}  is  obtained  for  a particular  frequency  (eigenvalue) 
of  oscillation,  the  velocity  components  (Equation  5)  may  be  obtained 
from  the  slope  of  £ , which  when  expressed  at  the  element  level  is 


3£ 

3x 

3£ 

3y 


G(Ue 


where  0 is  the  element  slope  matrix.  G is  defined  as: 


(it) 


G 


3N.  3N.  3N 

i J.  k 

3x  3x  3x 


3N.  3N . 3N, 
i_  ■)  k 

3y  3y  3y 


(15) 


g 

Thus,  { f, } having  been  calculated  and  the  shape  function  1 N ] known, 
the  horizontal  velocity  at  any  point  (x,y)  within  an  element  can  be 
determined  from  Equations  it  and  15.  For  this  study,  the  velocity  was 
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evaluated  at  the  centroid  of  each  element. 

15.  The  numerical  calculations  required  for  the  finite  element 
procedure  were  performed  on  a 0-635  Honeywell  computer  with  plotting 
performed  on  a Oalcomp  drum  plotter.  The  inputs  required  by  the  finite 
element  code  are  as  follow: 

a_.  The  total  number  of  elements  and  nodes, 
b^.  The  arrangement  of  nodes  for  each  element. 

£.  The  x,y  coordinates  of  each  node. 
d_.  The  depth  at  each  node. 
e_.  The  boundary  conditions. 

_f . The  digitized  natural  boundary  of  the  lake. 

The  basic  outputs  of  the  finite  element  code  are  as  follow: 

a_.  M natural  frequencies  (or  eigenvalues)  of  the  lake. 

b.  Normalized  surface  displacements  at  every  node  and 
centroid  of  each  element  for  each  natural  frequency. 

c_.  Normalized  x and  y velocity  components  at  the  centroid 
of  each  element  for  each  natural  frequency. 

Thus,  results  foi’  the  local  region  of  interest  (Figure  2)  consist  of 
normalized  surface  elevations  at  each  point  labeled  with  an  alpha  N 
or  alpha  E and  normalized  velocity  components  at  each  point  labeled 
with  an  alpha  E . 

Discretization  of  Lake  Erie 

16.  The  surface  area  of  Lake  Erie  (Figure  3)  was  approximated  by 
238  elements  with  264  total  nodal  points.  It  should  be  noted  that  many 

of  the  elements  are  quadrilateral  in  shape,  and  these  elements  were 

* 

represented  in  the  finite  element  code  by  four  CST's.  A local  region 
of  interest  near  Cleveland,  Ohio,  is  identified  by  a heavy  lined  bound- 
ary in  Figure  3.  This  area  is  expanded  as  Figure  2 for  element  defi- 
nition of  shape  with  nodes  and  centroids  for  elements  of  particular 
interest  being  assigned  numbers  for  future  identification  in  interpret- 
ing the  results.  The  choice  of  the  element  shapes  in  the  vicinity  of 


* J.  F.  Abel,  Memorandum  for  Record,  U.  S.  Army  Engineer  Waterways 
Experiment  Station,  CE,  Vicksburg,  Miss. 
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Figure  3»  Finite  element  representation  of  Lake  Erie 
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the  proposed  jetport  site  was  dictated  by  the  configurations  of  the  off- 
shore jetport  island,  its  shoreward  extensi  n,  and  s.hore  connection. 
Because  their  location  and  conf iguration  had  not  been  selected  at  the 

time  of  this  study,  tentative  configurations  for  this  study  were  chosen 

# 

based  on  available  prefeasibility  information.  Figure  h is  a schematic 
of  these  tentative  configurations.  The  finite  < Lement  representation  of 
these  assumed  configurations  and  locations  of  the  jetport  island  are 
shown  in  Figure  5a,  the  shoreward  extension  in  Figure  5b,  and  two  shore 
connections  in  Figures  5c  and  d.  It  is  emphasized  that  these  four 
configurations  and  locations  are  tentative  for  this  study  and  are  sub- 
ject to  change  as  LERTA's  feasibility  study  continues. 


Figure  A.  Schematic  of  tentative  configurations 


17.  The  depth  of  the  lake  is  represented  by  discrete  values  at 
each  nodal  point.  Figure  6 depicts  the  representative  depths  at  the 
nodal  points.  These  depths  were  taken  from  the  U.  S.  Lake  Survey 

* "The  Lake  Erie  International  Jetport  Project,"  Pre-Feasibility  Tech- 
nical Report,  March  1971,  Greater  Cleveland  Growth  Association. 
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Finite  element  representation  for  tentative  configurations 


Chart  L.S.j  which  is  referenced  to  a low  water  datum  of  570.5  ft.* 

18.  The  volume,  surface  area,  and  mean  depth  for  the  discretized 
representation  of  Lake  Erie  are  represented  in  the  following  tabulation. 


Parameter 

Volume,  cu  ft 
Surface  area,  sq  ft 
Mean  depth,  ft 


Observed 

16.8  * io12 

27.6  x 1010 

60.7 


Computed 

16.1*3  x 1012 
27.07  x io10 
60.69 


For  comparison,  the  factual  (observed)  values  of  these  parameters  for 
Lake  Erie  are  also  listed.  The  computed  values,  in  particular  the  mean 
depth,  agree  very  well  with  the  observed  values. 


Ver i f ication 


19.  To  verify  the  application  of  the  numerical  procedure  outlined 
in  the  previous  sections  to  the  investigation  of  the  effects  of  a jet- 
port  on  the  natural  oscillations  of  Lake  Erie,  the  natural  periods  and 
corresponding  normal  modes  of  variation  in  surface  elevations  for  the 
first  five  modes  of  oscillations  were  calculated  for  Lake  Erie  without 
the  jet port  structure.  The  calculated  natural  periods  of  the  free 
oscillations  are  shown  in  the  following  tabulation. 


Mode 


Computed  Observed 

Periods,  hr  Periods,  hr* 


1 

ll*  .1*3 

lU . 38 

2 

9.22 

9.ll» 

3 

6.01 

5.93 

1* 

1* . 31 

lt.15 

5 

3.87 

f From 

Reference  9. 

* A table  of  factors  for  converting  metric  (Si)  units  of  measurement 
to  U.  3.  customary  units  and  U.  S customary  units  to  metric  (Ol) 
units  is  given  on  page  4. 
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Also  tabulated  are  the  observed  periods  for  the  free  oscillations  of  Lake 

Erie  as  determined  from  spectra  of  records  of  lake  levels  by  Platzman  and 

9 

Rao.  The  discretized  depths  in  this  study,  discussed  in  para, -rush  17, 
were  corrected  to  correspond  to  actual  depths  of  the  lake  at  time  of 
observed  oscillations  by  considering  differences  in  the  reference  datum 
of  Chart  L.S.3  and  mean  lake  level  at  times  of  observation.  There  is 
relatively  good  agreement  between  the  observed  periods  and  the  periods 
computed  by  the  finite  element  procedure.  In  addition  to  the  comparison 
of  periods,  the  computed  amplitudes  of  the  relative  surface  variation 
for  the  fundamental  mode  of  seiching  for  13  locations  around  the  lake 
were  compared  with  the  observed  values  from  Platzman  and  Rao.10  The 
relative  amplitudes  around  the  lake  from  this  seiche  study  are  normal- 
ized amplitudes  such  that  actual  seiche  amplitudes  around  the  lake  can 
be  determined,  if  the  seiche  amplitude  at  one  location  is  known,  by 
multiplying  the  relative  amplitudes  by  the  known  seiche  amplitude.  For 
example,  actual  seiche  amplitudes  around  the  lake  for  a 3-ft  amplitude 
at  Toledo  can  be  estimated  by  multiplying  this  3-ft  amplitude  by  the 
relative  amplitudes  around  the  lake.  For  this  case,  the  actual  seiche 
amplitude  of  the  first  mode  at  Cleveland  would  be  3.0  times  0.27  or 
0.8  ft.  Table  1 shows  this  comparison,  which  exhibits  good  agreement 
around  the  lake.  The  maximum  deviation  between  computed  and  observed 
relative  amplitudes  occurs  near  Port  Clinton.  This  deviation  may 
reflect  the  reliability  of  the  observations  at  Port  Clinton  or,  more 
probable,  is  associated  with  the  numerical  representation  of  the  islands, 
topography,  arid  shoreline  in  the  eastern  basin  of  the  lake. 

Model  Results 


20.  For  the  full  lake  and  nearshore  region  without  a Jetport,  mode 
configurations  of  relative  amplitudes  and  depth-averaged  horizontal 
velocities  are  numerically  computed  for  the  first  five  modes  (para- 
graph 19)  or  eigenvalues.  Results  of  these  computations  for  the  full 
lake  are  presented  in  Figures  7-l6.  For  each  spatial  location  in  the 
nearshore  region,  the  depth-averaged  horizontal  velocity  is  defined  as 
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without  jetport  (normalized  to  maximum  elev: 
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jetport  (normalized  to  maximum 


a velocity  with  constant  direction  and  magnitude  over  the  erit ire  water 

column  and  whose  constant  direction  and  magnitude  represent  the  averaged 

velocity  over  the  entire  water  column.  In  addit  i n,  th  ma  jnil  ide  f 

the  depth-averaged  horizontal  velocity  is  a normalized  variable  like 

the  relative  seiche  amplitude  r<  fer  to 'paragraph  19).  The  velocity 

magnitudes,  which  are  computed  in  this  study,  are  normalized  to  a maxi- 
» 

mum  sei  :he  amplitude  of  1.0  ft.  For  the  first  mode  (fundamental),  the 
magnitude  of  the  normalized  velocity  vectors  in  the  Cleveland  area  must 
be  multiplied  by  the  actual  seiche  -unplitude  at  Toledo  to  obtain  the 
actual  velocity  magnitude.  However,  for  comparison  the  relative  magni- 
tudes can  be  used  without  scaling.  The  nearshore  results  without  a jet- 
port  for  the  first  five  modes  of  oscillation  are  presented  in  Figure  IT. 
For  four  generalized  configurations  and  locations  (Figure  5)  of  the  jet- 
port,  nearshore  results  of  the  relative  amplitudes  and  depth-averaged 
horizontal  velocities  for  the  first  five  modes  are  presented  in 
Figures  13-22.  The  calculated  periods  of  the  first  five  modes  of  free 
oscillations  for  different  jetport  configurations  are  presented  in 
Table  2. 


Conclusions  of  the  Seiche  Analysis 


21.  Based  upon  results  of  the  numerical  analysis  of  the  effects  of 
a jetport  offshore  near  Cleveland,  Ohio,  on  the  free  oscillations  of 
Lake  Erie,  several  ‘onclusions  are  apparent.  These  conclusions  hold 
for  all  five  modes  of  free  oscillation  unless  stated  otherwise  and  are 
as  follow: 

a.  Periods  of  the  first  five  modes  are  relatively  unchanged 
by  any  jetport  configuration.  There  exist  slight  in- 
creases (2. 5-5-0  min.)  in  the  periods  of  the  first  three 
modes  with  the  first  mode  exhibiting  the  largest  increase 
of  U-5  min. 

b^.  Changes  in  relative  seiche  amplitudes  and  horizontal 

velocity  are  confined  within  a distance  of  U-6  miles  from 
the  four  jetport  configurations. 
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Figure  18.  Mode  shape  and  normalized  current  for  first 
eigenvalue  rfith  different  .jetport  configurations 
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Figure  19 • Mode  shape  and  normalized  current  for  second 
eigenvalue  with  different  jetport  configurations 


Figure  20.  Mode  shape  and  normalized  current  for  third 
eigenvalue  with  different  jetport  configurations 
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Figure  21.  Mode  shape  and  normalized  current  for  fourth 
eigenvalue  with  different  jetport  configurations 
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c_.  Changes  iri  relative  seiche  amplitude  due  to  any  of  the 
four  jetport  configurations  are  equal  to  or  less  than 
0.1  ft  for  a 1.0- ft.  maximum  sci  'die  n:v,  1 i t.nd.  r :/  gr<*«»»  « r 
than  10  percent  change  in  seiche  amplitudes. 

(i.  Changes  in  relative  seiche  amplitude  and  horizontal 

gradients  of  this  amplitude  are  largest  for  jetport  con- 
figuration C and  smallest  for  jetport  configuration  A. 

e_.  The  largest  horizontal  velocities  are  generated  by  the 
fundamental  or  first  mode  of  oscillation,  which  is  the 
most  important  and  dominant  mode.  For  this  mode,  jetport 
configuration  A produces  a relative  velocity  of  0.5  fps 
(relative  to  1.0- ft  maximum  seiche  amplitude)  between 
shore  and  jetport  structure,  and  configuration  B produces 
a relative  velocity  of  1.0  fps  between  shore  and  jetport 
structure.  The  relative  velocity  in  this  area  without  a 
jetport  is  0.25  fps.  For  the  first  four  modes  of  oscil- 
lation, both  configurations  A and  B increase  the  hori- 
zontal velocity  between  the  jetport  structure  and  the 
shoreline  with  configuration  B producing  the  largest 
increases . 

_f.  Recurrence  intervals'*'  for  Buf falo-minus-Toledo  setup 

(seiche  wave  height)  have  been  computed  for  the  20-year 
period  19t0-1959.  Using  these  results,  the  expected 
frequency  of  occurrence  and  absolute  values  of  seiche 
velocity  for  the  fundamental  mode  were  estimated  for 
existing  conditions  and  configurations  A and  B (between 
the  jetport  structure  and  the  shoreline).  These  estimates 
are  presented  in  the  following  tabulation. 


Recurrence 

Seiche 

Velocity,  fps 

Interval,  yr 

Amplitude,  ft 

Ex i s t i ng 

Plan  A 

Plan  B 

0.25 

3 

0.75 

1.5 

3.0 

0.50 

3.5 

0.875 

1.75 

3.5 

1.0 

t 

1.0 

2.0 

t.O 

2 

5 

1.25 

2.5 

5.0 

5 

6 

1.5 

3.0 

6.0 

20 

7 

1.75 

3.5 

7.0 

g_.  Jetport  configurations  C ant.  D block  the  alongshore  flow 
near  the  Jetport  configurations  with  configuration  C 
producing  the  largest  decrease  in  these  horizontal 
velocities  and  creating  areas  near  the  shoreline  with 
minimum  or  no  circulation. 
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h_.  Of  the  four  generalized  Jetport  configurations  used  in  this 
study,  jetport  configuration  A has  the  least  effect  on 
seiche  period.-,  relative  seiche  amplitudes,  and  relative 
horizontal  velocities  for  the  first  five  modes  of  free 
oscillation  in  Lake  Erie. 

22.  This  numerical  study  was  conducted  to  evaluate  quantitatively 
the  effects  of  a jetport  in  Lake  Erie  near  Cleveland,  Ohio,  on  seiches 
(free  oscillations)  in  the  lake.  The  four  jetport  conf igurations  used 
in  this  study  are  tentative  and  general  in  configuration  and  location; 
however,  these  configurations  should  provide  results  from  which  the 
effects  of  a specific  configuration  can  be  extrapolated.  The  conclu- 
sions of  this  study  indicate  that  the  nearshore  effects  of  all  jetport 
configurations  are  minimal  for  seiche  periods  and  mode  configurations. 
However,  the  horizontal  velocity  regime  is  most  affected  by  the  Jetport 
configurations  with  the  island  configuration  (configuration  A)  having 
the  least  effect  on  the  velocity  regime.  These  results  for  the  velocity 
regime  indicate  that  additional  studies  are  needed  to  determine  the 
nearshore  effects  of  a jetport  structure  on  the  velocity  regime  for 
steady-state  and  time-dependent  wind-driven  circulation  in  Lake  Erie. 


4l 


PART  III:  TWO-DIMENSIONAL  DEPTH-AVERAGED 

HYDRODYNAMIC  MODELS 


General  Formulation 


I 


23.  The  generalized  three-dimensional  differential  equations 
governing  the  motion  of  an  infinitesimal  fluid  element  are  readily 
derived  from  basic  considerations  of  mass  and  momentum.  If  the  coordi- 
nate system  shown  in  Figure  23  is  used,  the  three-dimensional  governing 
equations  for  an  infinitesimal  fluid  element  can  be  expressed  in  the 
form  of  momentum  and  continuity  equations : 

Momentum  equations 
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Continuity  equation 
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w = velocity  in  z direction 
f = Coriolis  parameter 
p = density 


p = pressure 

K = coefficient  of  eddy  viscosity  or  diffusivity 
In  the  momentum  equations  the  first  term  is  the  local  acceleration 
term,  the  next  three  terms  are  the  convective  or  advective  acceleration 
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terms,  the  next  term  is  the  Coriolis  force  term,  followed  by  the 
pressure  gradient  term  and  the  viscous  stress  or  diffusion  terms.  The 
continuity  equation  is  simply  a statement  of  the  conservation  of  mass. 

2U.  Various  simplifying  conditions  are  applied  to  the  three- 
dimensional  equations  to  yield  a two-dimensional  system  of  equations 
that  can  be  more  readily  treated  by  numerical  techniques.  If  applica- 
tion is  made  only  to  problems  consistent  with  restrictions  imposed  by 
the  various  simplifying  assumptions,  then  valuable  information  can  be 
derived  from  two-dimensional  considerations. 

25.  A general  procedure  for  reducing  the  three-dimensional 
equations  to  a two-dimensional,  depth-averaged  formulation  will  be 
presented.  Most  two-dimensional,  depth-averaged  models  are  very  similar 
with  the  differences  due  to  the  specific  formulation  used  for  the  fric- 
tion and  diffusion  terms. 

26.  It  is  normal  to  simplify  the  continuity  equation  immediately 
when  working  with  "basically  incompressible  fluids."  The  Boussinesq 
approximation  is  made,  which  assumes  that  density  variations  are  small 
and  can  be  neglected  except  in  the  gravity  term  of  the  momentum  equation. 
Thus,  even  if  the  fluid  is  not  considered  as  strictly  incompressible, 
the  variation  in  density  is  considered  as  small  and  the  incompressible 
form  of  the  continuity  equation  is  used.  Two-dimensional  formulations 
almost  always  consider  the  flow  to  be  well  mixed  and  the  density  to  be 
constant.  Equation  19  reduces  to: 
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27.  The  vertical  component  of  the  local  acceleration,  advective 
acceleration,  and  velocity  is  assumed  to  have  a negligible  effect  on 
the  momentum  equations.  This  reduces  the  momentum  equation  to: 
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Equation  23  is  now  simply  a statement  of  the  hydrostatic  variation  of 
pressure  if  density  is  a constant.  The  hydrostatic  pressure  assumption 
is  reasonable  for  shallow-water  waves,  i.e.  long-period  waves,  where  the 
pressure  fluctuations  due  to  a wave  are  transmitted  virtually  unattenu- 
ated from  the  surface  to  the  bottom.  For  intermediate-  or  short-period 
waves  the  hydrostatic  assumption  may  not  be  a good  representation  of 
the  pressure. 

28.  Equations  20-23  are  now  in  the  appropriate  forms  for  applying 
the  depth  averaging  process.  Integrating  the  continuity  equation. 
Equation  20,  in  the  z direction  over  the  water  depth  gives 
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At  the  water  surface. 


z = n(x,y,t) 


must  be  a streamline.  In  other  words. 


(2U) 


dt 
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which  becomes 
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Therefore,  the  free-surface  boundary  condition  is 
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If  the  so-called  rigid  lid  assumption  is  made,  then  the  boundary 
condition  at  the  surface  reduces  to 


w 5 0 


(26) 


The  boundary  condition  at  the  bottom,  z = -h  , is  similar  to  that  at 
the  surface;  however,  since 
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the  boundary  condition  reduces  to 
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Substituting  Equations  25  and  27  into  2U  yields 
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Using  differentiation  of  an  integral  involving  a parameter  in  the  limits 
of  the  integral. 


U6 


9_ 

3x 


n 


n(x,y) 


u dz  = u 


9n 

9x 


+ u 


9h 


9u 


+ T-  dz 


9x 


9x 


-h(x,y) 


z=n 


z=-h 


-h 


or 


r\u 

— dz  = 
9x 

9 

3x 

ru 

1 u dz  - u 

9h 

- u 
9x 

Lv, 

J 

1 

z=n 

9h 

9x 


z=-h 


(29) 


In  a similar  manner 


(30) 


Substituting  Equations  29  and  30  into  28  and  cancelling  terms  gives 


n n 


The  u and  v components  of  velocity  are  assumed  to  be  of  the  form 


u(x,y,z,t)  = U(x,y ,t ) [ 1 + u'(z)] 


v(x,y,z,t)  = V(x,y,t)[l  + v'(z)] 


where 


dz  = 0 


dz  = 0 
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U = x-component  of  depth-averaged  horizontal  velocity 
u'  = depth-dependent  perturbation  in  horizontal  velocity  in  x 
direction 

V = y-component  of  depth-averaged,  horizontal  velocity 
v'  = depth-dependent  perturbation  in  horizontal  velocity  in  y 
direction 

If  U(x,y,t)  and  V(x,y,t  ) are  depth-averaged  velocities  given  by 


and 
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then  the  continuity  equation  reduces  to 
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29.  From  the  simplified  momentum  equation  in  the  z direction. 
Equation  23, 
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Integrating  from  the  water  surface  to  an  arbitrary  depth  z gives 


pressure  with'  horizontal  distance  at  the  free  surface.  These  terms  are 
normally  neglected,  and  the  pressure  gradient  is  given  by 
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30.  The  momentum  equation  in  the  x direction  can  now  be  inte- 
grated over  the  depth  after  substituting  from  Equation  36  for  the 
pressure  gradient. 
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Considering  the  individual  terms  in  this  equation  and  performing  the 
indicated  integration  gives 
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(tx)  dZ 


-h 


where  is  stress  in  the  x direction.  Therefore,  the  basic  momen- 

tum equation  in  the  x direction  becomes 
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where 

t = the  x-component  of  wind  stress  at  the  water  surface 

WX 

= the  x-component  of  bottom  friction 

x5X 

31.  In  a similar  manner  the  momentum  equation  in  the  y direc- 
tion becomes 
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where 


t = the  y-component  of  wind  stress  at  the  water  surface 
wy 
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t,  = the  y-component  of  bottom  friction. 

By 

32.  Equations  31* , 39  > and  ^0  form  the  basis  for  most  two- 
dimensional,  depth-averaged  models.  These  equations  can  be  formulated 
«in;termr  of  finite  difference  equations  for  solution.  It  should  be 
noted  that  while  the  equations  are  formulated  in  terms  of  two- 
dimensional  horizontal  velocity  components,  a pseudo-three-dimensional 
effect  is  present  since  actual  depths  are  input  to  the  model. 


Discussion  of  Specific  Models 

33.  The  two-dimensional  depth-averaged  formulation  was  shown  to 

reduce  to  the  system  of  equations  represented  by  Equations  3^,  39,  and 

ItO.  These  differential  equations  must  be  discretized  and  formulated  in 

terms  of  finite  difference  equations  for  solution  by  numerical  tech- 
12 

niques.  The  finite  difference  representation  of  the  differential 
equations  is  not  unique;  it  can  be  formulated  in  terms  of  backward 
differences,  forward  differences,  central  differences,  or  some  combi- 
nation of  these  difference  schemes.  The  problem  also  can  be  formulated 
such  that  it  is  explicit,  implicit,  or  implicit-explicit  in  time.  Each 
method  of  formulation  does  offer  advantages  and  disadvantages,  although 
some  formulation  schemes  have  been  shown  to  be  completely  unsatisfactory. 
In  general,  the  explicit  techniques  are  simpler  to  formulate  and  require 
less  computer  time  for  a given  size  time  step.  The  explicit  schemes, 
however,  have  stability  criteria  that  normally  require  a smaller  step 
size  than  the  comparable  calculation  using  an  implicit  scheme. 

3^.  Basically  the  two-dimensional  depth-averaged  models  are  used 
for  two  different  types  of  problems.  In  the  first  case  the  determi- 
nation of  circulation  patterns  is  the  primary  objective  with  wind  or 
tides  being  the  predominant  forcing  function.  The  second  type  of 
application  is  used  mainly  to  calculate  the  displacement  of  the  free 
surface  during  storm  surge.  Experience  has  shown  that  the  importance 
of  the  various  terms  in  the  governing  equations  and  the  boundary  condi- 
tions vary  with  the  application.  For  example,  the  nonlinear  advective 
terms  and  horizontal  diffusion  terms  are  normally  neglected  in  the 
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storm  surge  models  while  they  may  be  very  important  in  circulation 
models,  particularly  near  the  shore. 

35-  Several  specific  two-dimensional  depth-averaged  models  will 
now  be  discussed.  The  departure  of  each  model  from  the  basic  formu- 
lation of  the  problem  will  be  outlined.  Applications  and  verification 
of  the  model  will  be  discussed  along  with  the  general  conclusions  of 
the  authors.  The  models  discussed  are  representative  of  the  many  two- 
dimensional  models  available. 

Model  of  T.  J.  Simons 

1 3 

36.  The  two-dimensional  numerical  model  by  T.  J.  Simons  is  a 
straightforward  application  of  the  two-dimensional,  depth-averaged 
equations.  The  model  is  formulated,  however,  in  terms  of  volumetric 
flow  rate  per  unit  width  (foot)  rather  than  average  velocities. 
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(P,Q)  = J'  ( u , v ) dz 

-h 


(1*1) 


where  P and  Q are  volumetric  flows  in  x and  y directions, 
respectively.  The  depth-averaged  velocities  can  be  obtained  by  averag- 
ing the  volume  transport  over  depth. 
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37*  The  equations  used  in  the  model  to  calculate  volume  transport 
are  as  follows: 
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These  equations  are  completely  consistent  with  Equations  3U,  39,  and  1*0. 

38.  Simons  expresses  the  bottom  stress  in  terms  of  a bottom 
friction  coefficient  and  the  components  of  volume  transport. 
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The  bottom  stress  coefficient  is  considered  to  be  represented  by  the 
following  simple  relationships: 


a_. 

B = a/h  where 

a = 0.01 

cm/ sec 

b. 

2 

B = b/h  where 

b ~ 100 

2, 

cm  /sec 

c. 

B = c|P,Q|/h2 

where  c 

= 0.0025 

39-  The  external  forces  during  circulation  are  the  surface  pres- 
sure and  the  wind  stress.  In  these  calculations  the  surface  pressure 
is  set  equal  to  zero.  The  surface  wind  stress  is  related  to  the  wind 
velocity;  however,  this  relationship  is  not  considered  in  this  model. 
Instead,  more  or  less  typical  wind  stress  fields  are  specified.  The 
wind  stress  fields  are  taken  to  simulate  the  general  behavior  of  the 
wind  following  the  passage  of  an  atmospheric  front.  The  winds  tend  to 
increase  sharply  for  a relatively  short  period  of  time  and  thereafter 
vary  more  slowly.  At  the  same  time  the  wind  field  will  move  across  the 
lake  in  the  general  direction  of  the  wind  stress  itself.  Idealizing 
this  situation,  a semi-infinite  stress  band  moving  with  a constant 
translation  speed  and  having  a linear  increase  of  intensity  of 

wind  stress  over  a period  of  time  T^  and  a constant  intensity  after 

time  T is  used  in  the  model, 
s 

1*0.  The  system  of  equations  is  solved  using  central  differences 
and  an  implicit  scheme.  A space-staggered  lattice  is  used  on  which 
velocity,  water  level  displacement,  and  water  depth  are  described  at 
different  locations  within  a grid  cell. 
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Ul.  The  primary  application  of  the  model  is  to  predict  circulation 
patterns  in  lakes.  The  model  has  been  applied  to  Lake  Ontario  both  with 
and  without  the  nonlinear  advective  terms  in  the  equation.  A 5-km 
rectangular  grid  was  used,  and  circulation  patterns  are  shown  at  various 
times  during  the  passage  of  a hypothetical  wind  field. 

U2.  Based  upon  numerical  experimentation  with  the  various  param- 
eters in  the  model,  the  following  conclusions  as  to  the  effects  of  these 
parameters  are  presented: 

a..  The  actual  effect  of  the  nonlinear  terms  does  not  seem  to 
justify  the  tremendous  increase  in  computational  effort. 

b_.  The  nonlinear  effects  are  completely  obscured  by  the 
effects  of  bottom  topography,  bottom  friction,  and 
horizontal  diffusion. 

c_.  Increased  efforts  should  be  directed  at  increasing  hori- 
zontal resolution  and  improving  numerical  treatment  of  the 
boundari es . 

d_.  The  actual  formulation  used  for  the  bottom  stress  coef- 
ficient does  not  have  a great  effect  on  the  circulation 
pattern  as  long  as  a reasonable  formulation  is  used. 

e_.  The  smoothing  effect  of  the  horizontal  diffusion  term 
is  considerable. 

_f . The  vertically  integrated  or  averaged  velocities  produced 
by  the  model  are  not  indicative  of  the  actual  velocities 
found  in  a lake.  In  general,  the  volume  transport  model 
tends  to  underestimate  the  magnitude  of  the.  surface 
velocities,  the  time  variation  of  the  flow,  and  the 
Coriolis  force  effect. 

£.  A three-dimensional  or  layered  model  is  needed  for 
velocity  distributions  in  a lake. 

Model  of  George  W.  Platzman 

i lU 

43.  The  primary  application  of  the  model  by  Platzman  is  to 

predict  surface  displacement  during  storm  surge  on  lakes.  The  formu- 
lation of  the  problem  is  in  term  of  volume  transport  rather  than  average 
velocity.  The  basic  equations  are  essentially  identical  with  those  by 
Simons,  Equations  U3-^5,  except  that  the  nonlinear  advective  terms  and 
the  horizontal  diffusion  terms  are  neglected. 

Ub.  While  the  equations  used  by  Platzman  are  similar  to  those  used 
by  Simons , tne  wind  stress  and  bottom  friction  terms  are  handled 


55 


differently.  Coefficients  that  are  functions  of  the  Ekman  number  are 

introduced  in  such  a manner  that  the  solution  for  steady-state  transport 

15 

reduces  to  the  analytical  solution  obtained  by  Ekman  for  surface  slope 
transport  and  wind  transport.  Basically,  the  bottom  friction  is  related 
to  the  transport,  and  wind  stress  is  related  to  the  square  of  the  wind 
speed;  however,  the  various  parameters  are  expressed  as  a function  of 
the  Ekman  number. 

1*5-  The  system  of  equations  is  solved  using  central  differences 
in  time  and  space  and  a space-staggered  lattice  upon  which  velocities, 
water  level  displacements,  and  water  depths  are  described  at  different 
locations  within  a grid  cell.  The  model  provides  for  a time  and  spatial 
variation  of  wind  stress  by  input  of  actual  wind  data. 

46.  The  model  is  applied  to  nine  actual  extreme  storm  surge  con- 
ditions on  Lake  Erie  with  wind  stress  information  based  upon  hourly 
surface  wind  observations  at  six  weather  stations  located  on  the 
periphery  of  the  lake.  The  wind  field  on  the  lake  varies  with  time  and 
position  and  is  based  upon  a weighted  average  of  the  wind  data  from  the 
six  weather  stations.  Computed  lake-level  configurations  have  been 
compared  with  observed  data  for  the  nine  extreme  storm  surge  conditions 
considered. 

47.  The  following  conclusions  on  the  applicability  of  the  model 
for  storm  surge  analysis  are  presented: 

a^.  The  computed  wind  setup  values  agree  well  with  observed 
values . 

b.  Perhaps  the  most  conspicuous  defect  of  the  computations 
is  failure  to  predict  the  free  oscillations  of  the  lake 
after  the  passing  of  the  storm. 

c_.  Improvements  in  the  model  would  probably  result  from  more 
careful  consideration  of  the  effect  of  differences  in 
anemometer  exposure  for  the  various  wind  stations  used  in 
the  analysis. 

Model  of  J,  J.  Leendertse 

48.  The  primary  application  of  the  model  by  Leendertse'*’^  is  for 
circulation  patterns  for  tidal  flows  in  bays  and  estuaries.  The  basic 
governing  equations  are  identical  with  Equations  34,  39*  and  40.  The 
horizontal  diffusion  terms  are  neglected;  while  wind  stress  is  included 
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in  the  formulation,  tidal  elevations  are  the  primary  forcing  function. 
The  bottom  stress  term  is  considered  proportional  to  the  squared  veloci- 
ty and  is  formulated  in  terms  of  the  Chezy  coefficient,  C,  as  follows: 


(US) 


(U9) 


The  wind  stress  is  formulated  in  the  normal  manner  in  terms  of  the 
square  of  the  wind  speed. 

U9.  An  explicit-implicit  scheme  is  used  to  solve  the  system  of 
equations  using  a space-staggered  grid  arrangement  where  velocities, 
surface  displacement,  and  depths  are  defined  at  different  locations  in 
the  grid  system.  This  model  has  been  applied  to  the  North  Sea,  Jamaica 
Bay,  Tokyo  Bay  and  other  large  bodies  of  water.  The  general  conclusions 
regarding  the  model  are  as  follows : 

si.  For  applications  to  bays,  harbors,  inlets,  etc.,  where 
tidal  elevations  are  the  predominant  forcing  function 
and  where  a significant  volume  transport  occurs,  the 
depth-averaged  circulation  patterns  appear  to  give  reason- 
able results. 

b_.  The  model  can  also  be  used  to  model  water  waves  generated 
by  nuclear  explosions,  tsunami  waves,  and  other  long- 
period  wave  motion. 

50.  This  model  has  no  direct  applicability  to  the  problem  of  wind- 
driven  circulation  or  storm  surge  on  Lake  Erie.  The  implicit-explicit 
numerical  scheme  used  to  solve  the  system  of  difference  equations  was 

of  primary  interest  in  this  model.  This  basic  numerical  procedure  has 
several  desirable  features  that  might  prove  beneficial  if  the  procedure 
were  applied  to  a model  of  Lake  Erie.  This  particular  solution  scheme 
was  investigated  since  it  had  not  been  used  previously  in  a lake  model. 
Model  of  Frank  D.  Masch 

IT 

51.  The  model  by  Masch  is  primarily  intended  for  application  to 
estuaries  where  significant  tidal  flats  exist.  The  basic  equations  are 
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formulated  in  terms  of  volume  transport  and  are  essentially  identical 
to  Equations  it 3-^*5*  The  bottom  friction  term  is  expressed  in  terms  of 
the  Manning  friction  factor,  n : 
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Horizontal  diffusion  is  neglected.  This  model  is  basically  an  extension 

l8 

of  the  Reid-Bodine  model. 

52.  The  model  contains  provisions  for  flooding  of  tidal  flats, 
such  obstructions  to  flow  as  weirs  or  breakwaters,  and  other  phenomena 
somewhat  peculiar  to  tidal  flows  in  estuaries.  Several  computational 
boundaries  can  be  represented  in  the  model  indicating  ocean  boundaries, 
river  inflows,  or  artificial  boundaries  established  in  the  estuary  for 
computational  convenience.  This  model  has  been  applied  to  several 
locations  including  Matagorda  Bay  and  Galveston  Bay,  Texas,  and 
Masonboro  Inlet,  North  Carolina.  As  a result  of  these  applications,  it 
was  concluded  that: 

a.  Surface  displacements  and  volume  transport  due  to  tidal 
forces  can  be  predicted  reasonably  well  with  the  model. 

b.  Flooding  of  tidal  flats,  underwater  obstructions,  break- 
waters, and  similar  obstruction  to  flow  can  be  included 
in  the  model  if  step  sizes  are  small  enough  for  the 
necessary  spatial  resolution. 

53-  This  model  has  no  direct  application  to  the  problem  of  wind- 
driven  circulation  or  storm  surge  on  Lake  Erie.  This  model  was  investi- 
gated because  of  the  detailed  treatment  of  numerical  procedures  for 
treating  breakwaters,  weirs,  underwater  obstructions,  flows  from  streams 
into  the  region  being  considered,  flooding  regions,  and  other  hydro- 
dynamic  features.  These  phenomena  have  not  been  treated  in  detail  in 
previous  lake  models,  and  their  use  may  be  necessary  in  modeling  a 
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jetport  in  Lake  Erie  and  its  effect  on  circulation,  particularly  in 
Cleveland  Harbor. 

Model  of  David  F.  Paskausky 

5^.  The  model  by  Paskausky'*'*  has  its  origin  in  the  three- 
dimensional  equations  of  hydrodynamics;  however,  the  development  of  the 
model  differs  significantly  from  the  models  discussed  previously.  The 
primary  application  of  the  model  is  for  storm  surge;  however,  the 
volume  transport  per  unit  width  (foot)  is  the  model  output.  A rigid 
lid  assumption  is  made  at  the  water  surface. 

55*  The  development  of  this  model  initially  follows  the  basic 
formulation  of  the  two-dimensional  depth-averaged  model.  Consider 
Equations  20-23,  which  are  the  simplified  three-dimensional  differential 
equations  prior  to  vertical  integration  over  the  water  depth.  Differ- 
entiation of  Equation  21  by  - 9/3y  and  Equation  22  by  9/9x  introduction 

of  the  vorticity  component  £ and  shear  stress  components  t and  t 

x y 


9v  9u 
9x  3y  ’ Tx 


and  addition  of  the  two  equations  yield  the  vorticity  equation: 


(52) 


In  Equation  52,  A is 
advection  term,  C is 
advection  term,  E is 
torque  per  unit  volume 


the  topographic  term,  B is  the  horizontal 
the  planetary  vorticity  term,  D is  the  vertical 
the  lateral  friction  term,  and  F is  the  net 
exerted  on  the  top  and  bottom  of  a water  parcel. 
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The  vertical  advection  term  vanishes  for  a well  mixed,  constant-density 
application . 

56.  By  combining  terms  and  ignoring  the  vertical  advection  term. 
Equation  52  can  be  written  in  the  form: 
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where 

Ah  = i (a/9x)  + j (a/3y) 
vh  = iu  + Jv 
v2  = (a2/ax2)  + (a2/ay2) 
curl  t = (8t  /3x)  - (3t  /3y) 

y ^ 


in  which  i is  unit  vector  in  the  x direction  and  j is  unit  vector 
in  the  y direction.  Equation  53  is  now  integrated  over  the  depth,  D , 
to  obtain  a two-dimensional  depth-averaged  model. 
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Considering  the  individual  terms  in  the  equation  for  the  well  mixed, 
constant-density  condition  and  assuming  the  bottom  friction  stress 
equal  to  the  mass  flux  multiplied  by  a bottom  friction  coefficient  o 
gives 
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Where  is  the  horizontal  diffusion  coefficient.  Equation  51*  reduced  to 


+ Vu  • (f  + C)  V = K,  V‘.C  - a?  - ^ curl  t 


3t  h 


h “h  h- 


(55) 


An  Ekman  boundary  layer  approximation  for  the  bottom  friction  coeffi- 
cient is  used  at  the  sea  bed: 


0.5K  f\1/2 

s3-) 


where  K is  the  vertical  diffusion  coefficient, 
v 

57.  Integrating  the  continuity  equation  (Equation  1*5)  across  the 


depth  yields 

D 


o 


Using  the  "rigid  lid"  boundary  condition  (Equation  26)  at  the  free 
surface  and  the  bottom  boundary  condition  (Equation  27)  reduces  to 


sir (Du)  * ■£  {In) - 0 


(56) 


A volume  transport  stream  function  ip  exists,  such  that 
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In  terms  of  <|>  the  vorticity  component  £ is  given  by 
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(57) 


(58) 
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Use  of  Equation  58  with  a knowledge  of  C,  from  Equation  55  allows 
evaluation  of  ip  by  relaxation  subject  to  appropriate  boundary  condi- 
tions. The  velocity  field  can  then  be  calculated  from  Equation  57. 

58.  The  model  has  been  applied  to  Lake  Ontario  for  a simplified 

storm  condition  passing  along  the  long  axis  of  the  lake.  The  velocity 

conditions  in  the  lake  are  calculated  at  various  times  during  the 

passage  of  the  storm.  The  model  has  also  been  applied  to  Lake  Erie 

20 

during  ttie  passage  of  storm  Agnes. 

59.  While  the  results  from  this  model  are  qualitatively  interest- 
ing, certain  of  its  assumptions  are  subject  to  question.  The  volume 
transport,  which  is  integrated  over  the  depth  of  the  lake,  and  hence 
the  surface  elevations  may  be  reasonably  correct;  however,  these  values 
cannot  be  related  to  prototype  velocity  values  in  the  lake,  which  are 
highly  depth  dependent.  The  applicability  of  the  rigid  lid  assumption 
made  in  this  model  is  questionable  for  storm  surge. 

Limitations  of  the  Models 


60.  From  the  literature  survey  of  the  various  two-dimensional 
hydrodynamic  models  available  and  based  upon  the  experience  at  WES  with 
some  of  these  models,  the  applicability  and  limitations  of  these  models 
become  fairly  apparent: 

si.  Two-dimensional  models  will  probably  yield  reasonable 

results  in  lakes  for  storm  surge  (free  surface  displace- 
ment) if  a free  surface  rather  than  a rigid  lid  boundary 
condition  is  used  at  the  lake  surface.  Free  surface 
displacement  is  basically  a function  of  the  total  mass 
transport,  and  two-dimensional  models  appear  reasonably 
capable  of  representing  this  phenomenon. 

b.  Two-dimensional  models  will  not  yield  results  for  depth- 
dependent,  wind-driven  lake  currents,  which  are  repre- 
sentative of  the  currents  existing  in  the  lake. 

c^.  Two-dimensional  models  probably  yield  reasonable  values 
for  both  surface  displacements  and  currents  in  estuaries 
and  bays  if  large  salinity  variations  are  not  present. 

61.  The  first  and  second  conclusions  appear  in  conflict  at  first 
inspection.  This  apparent  conflict  can  be  resolved  by  considering  the 
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difference  in  surface  elevations  required  to  produce  a specified  veloc- 
ity. Very  small  changer  in  surface  elevation  can  produce  significant 
velocities.  Thun,  the  two-dimensional  mathematical  model  can  predict 
values  of  surface  displacement  in  reasonable  agreement  with  experimental 
data  if  the  total  mass  transport  is  approximately  correct.  At  the  same 
time  very  poor  representations  of  the  actual  currents  in  the  lake  are 
produced  by  the  models  due  to  the  complex  depth-dependent  velocities 
existing  in  the  lake.  In  most  freshwater  lakes  the  throughflow  is 
small,  and  wind  is  the  major  current-producing  forcing  function.  The 
current  pattern  is  extremely  depth  dependent,  varying  in  direction  from 
surface  to  bottom  by  essentially  180  deg.  A depth-integrated  velocity 
such  as  produced  by  the  two-dimensional  models  will  generally  yield  a 
very  small  average  velocity  proportional  to  the  lake  throughflow  and 
completely  fail  to  indicate  the  very  complex  depth-dependent  velocities 
existing  in  the  lake. 

62.  Since  the  two-dimensional  models  have  no  provisions  for  a 
return  flow  along  the  bottom  of  the  lake,  they  must  establish  return 
flow  regions  in  the  lake  itself  when  wind  is  the  primary  forcing  func- 
tion. The  depth-averaged  flows  predicted  by  the  two-dimensional  models 
will  tend  to  be  in  the  wind  direction  in  the  shallow  areas  of  the  lake 
with  a return  flow  in  the  deeper  portions.  This  flow  further  distorts 
the  current  pattern  obtained  from  the  two-dimensional  model  from  those 
actually  existing  in  the  lake. 

63.  Because  conditions  in  estuaries  and  bays  normally  are  more 
consistent  with  tVie  basic  assumptions  used  in  the  derivation  of  two- 
dimensional  models,  the  results  are  in  better  agreement  with  prototype 
conditions.  For  estuaries  and  bays,  the  primary  forcing  function  is 
the  tide,  and  significant  mass  transport  does  occur.  The  velocity 
distribution  over  the  depth  is  generally  in  approximately  the  same 
direction  and  of  the  same  order  of  magnitude.  A depth-averaged  veloci- 
ty therefore  has  some  physical  meaning  in  the  estuary  or  bay,  and  nu- 
merical velocity  results  compare  more  favorably  than  was  the  case  in 
lakes . 
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PART  IV:  THREE-DIMENSIONAL  HYDRODYNAMIC  MODELS 


General  Formulation 


6h.  The  three-dimensional  hydrodynamic  models  have  only  recently 
been  developed  to  the  state  where  application  to  complex  geometries  with 
reasonable  resolution  appears  possible.  The  delay  in  developing  three- 
dimensional  models  can  be  attributed  to  the  unavailability  of  extremely 
large,  fast  computers  as  well  as  to  formulation  difficulties.  As  previ- 
ously indicated,  the  depth-averaged  two-dimensional  models  have  been 
extensively  applied  to  the  hydrodynamics  of  lakes  and  oceans.  The  two- 
dimensional  models  have  provided  valuable  quantitative  information  for 
certain  types  of  problems;  however,  for  many  other  problems  the  two- 
dimensional  models  yield,  at  best,  qualitative  results.  Certainly  the 
many  pressing  problems  associated  with  pollutant  and  sediment  transport 
require  the  details  of  the  three-dimensional  velocity  profile  in  the 
lake,  and  such  information  cannot  be  obtained  from  the  depth-integrated 
two-dimensional  system  of  equations.  The  shortcomings  of  the  depth- 
averaged  models  are  most  apparent  when  the  bottom  shear  stress  is  con- 
sidered. In  the  two-dimensional  models  the  shear  stress  is  normally 
related  to  the  net  horizontal  mass  flux.  However,  for  major  portions 
of  a lake,  the  flow  at  the  bottom  may  be  a return  flow  that  is  not  in 
the  direction  of  the  net  mass  flux. 

65.  The  three-dimensional  models  can  be  broadly  classified  into 
the  following  categories: 

a, .  Constant-density,  steady-flow  models 

b.  Constant-density,  unsteady-flow  models 
c_.  Variable-density,  unsteady-flow  models 

Within  each  category  the  models  can  vary  considerably  depending  upon 
the  assumptions  that  are  imposed.  Some  of  the  assumptions  are  common 
to  most  of  the  models;  others  pertain  to  a specific  model.  A variable- 
density,  steady-flow  result  can  be  obtained  from  model  c_  and  is  not 
listed  separately  since  no  basically  different  formulation  has  been 
developed  for  problems  of  this  type.  The  constant-density,  steady-flow 
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models  can  be  formulated  more  simply  than  the  constant-density, 
unsteady-flow  models. 

66.  The  three-dimensional  governing  equations  for  an  infinitesimal 
fluid  element  have  been  given  previously  in  Equations  16-19-  The  basic 
equations  for  the  three-dimensional  models  are  obtained  from  these  equa- 
tions subject  to  certain  assumptions  that  are  common  to  most  models. 

67.  The  Boussinesq  approximation  is  normally  made.  This  approxi- 
mation assumes  that  density  variations  are  small  and  can  be  neglected 
except  in  the  gravity  term  of  the  momentum  equation.  Thus,  even  if  the 
fluid  is  not  strictly  treated  as  incompressible,  the  continuity  equation 
is  reduced  to  the  incompressible  form.  In  addition,  for  most  existing 
models,  the  vertical  component  of  the  local  and  advective  acceleration 

is  assumed  to  have  a negligible  effect  on  the  vertical  momentum  equation. 
This  assumption  reduces  the  vertical  momentum  equation  to  the  equation 
for  hydrostatic  pressure  variation  in  a compressible  fluid. 

68.  The  basic  equations  from  which  the  development  of  most  three- 
dimensional  models  proceed  are  therefore  given  by 
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In  addition  to  these  equations,  an  energy  equation,  a constitutive 
equation  for  the  density  and  diffusion  equations  for  substances  of 
interest,  may  be  required,  depending  upon  the  particular  problem  formu- 
lation. The  K used  in  the  momentum  equation  is  generally  a 
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combination  of  the  molecular  viscosity  and  the  eddy  viscosity.  In 
practice,  however,  the  molecular  viscosity  is  often  neglected  since  it 
is  much  smaller  than  the  eddy  viscosity  for  lake  applications. 

Discussion  of  Specific  Models 

69.  The  three-dimensional  models  vary  considerably  in  their  method 
of  development;  therefore,  no  general  formulation  of  the  model  governing 
equations,  such  as  was  possible  in  the  two-dimensional  equation,  will  be 
attempted.  Instead,  specific  models  illustrating  the  various  model 
categories  will  be  considered  in  sufficient  detail  to  indicate  the 
formulation  procedure  as  well  as  the  strengths  and  weaknesses  of  the 
individual  models.  In  general,  the  application  of  three-dimensional 
models  to  actual  lake  geometries  is  still  in  its  infancy.  Considerable 
work  is  still  required  in  this  area,  and  many  of  the  questions  concern- 
ing the  accuracy  of  the  various  models  can  be  resolved  only  when  proto- 
type data  of  the  required  quality  and  quantity  become  available. 

Three-dimensional,  steady- 

state  , constant-dens : by 

model  of  R.  T.  Gedney  and  W.  Lick 

21  22 

70.  This  model  ’ is  based  upon  the  shallow-lake  theory  de- 
veloj cd  by  Welander.  The  model  can  account  for  irregular  boundaries  and 
variable  bottom  topography;  it  is,  however,  restricted  to  constant 
density  and  steady  flow.  The  primary  application  of  the  model  is  for 
calculation  of  three-dimensional  steady-state  velocities  when  ^he  lake 
is  in  an  unstratified  condition.  A rigid  lid  assumption  is  imposed  at 
the  surface  of  the  lake  to  filter  out  gravity  waves. 

71.  The  basic  model  is  capable  of  describing  flows  in  a lake 
where  the  depth  is  comparable  to  or  less  than  the  Ekman  layer  (shallow 
lake)  as  well  as  in  a deep  lake.  Several  assumptions  result  in  simpli- 
fication of  the  governing  equations  when  they  are  applied  to  shallow- 
lake  conditions.  The  model  is  three-dimensional;  however,  the  computer 
time  requirements  are  reasonably  small.  This  small  requirement  is 
possible  since  the  finite  difference  numerical  portion  of  the  model  is 
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actually  two-dimensional  with  the  variation  in  flow  parameters  along 
the  vertical  or  third  dimension  being  obtained  from  a closed-form  ana- 
lytical solution. 

7 2.  The  basic  starting  equations  for  the  model  are  those  listed 
in  Equations  59-62.  These  equations  are,  however,  considerably  simpli- 
fied by  assumptions  consistent  with  the  intended  application  to  shallow 
lakes : 


a . The  lake  is  considered  to  be  homogeneous,  resulting  in 
the  density  being  considered  as  a constant  even  in  the 
gravity  term  in  the  vertical  momentum  equation.  This 
assumption  of  homogenity  in  Lake  Erie  is  found  to  be  valid 
approximately  from  about  the  middle  of  fall  to  the  be- 
ginning of  summer. 


b_.  Tiie  nonlinear  inertial  forces  are  assumed  to  be  small  in 
comparison  with  the  Coriolis  force  and  are  neglected  in 
the  momentum  equations . The  Rossby  number  is  a measure 
of  the  ratio  of  the  nonlinear  inertia  forces  to  the 
Coriolis  force  and  is  defined  by  R^  = U^/fL  , where 
is  a characteristic  velocity  in  the  lake,  f is  the 
Coriolis  parameter,  ar.d  L is  a characteristic  horizontal 
scale.  A calculation  using  values  characteristic  of  Lake 
Erie  shows  that  except  in  localized  regions  where  the 
horizontal  length  scale  is  smaller  than  50  km,  the  non- 
linear inertia  terms  will  be  relatively  unimportant. 

These  terms  reduce  the  governing  equation  to  a linear 
system;  however,  the  solutions  obtained  from  the  model 
near  islands  and  in  shallow  shore  regions  should  be 
applied  with  caution. 


The  effects  of  the  horizontal  internal  friction  of  the 
fluid  are  ignored,  which  eliminates  the  terms  involving 
the  horizontal  eddy  viscosity.  The  horizontal  Ekman  number 
E^  is  a measure  of  the  relative  importance  of  the  hori- 
zontal diffusion  terms  by  comparison  with  the  Coriolis 
force  and  is  defined  by  Ejj  = Ajj/fL^  where  Ah  is  the 
horizontal  eddy  diffusion  coefficient.  Typical  values 
for  Lake  Erie  indicate  that  Eh  is  small  and  can  be 
ignored  except  in  the  boundary  layer  near  the  shoreline. 

The  boundary  layer  near  the  shore  is  normally  smaller, 
however,  than  the  typical  grid  size  used  in  the  calcu- 
lations . 


d.  The  pressure  is  assumed  to  be  hydrostatic  as  a result  of 
neglecting  vertical  acceleration  and  velocity  terms. 

This  standard  assumption  for  most  numerical  models  is 
probably  justified  except  near  the  shoreline  where  sig- 
nificant vertical  velocities  and  accelerations  can  be 
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present  as  a result  of  bringing  to  rest  the  velocity 
component  normal  to  the  shore. 

£.  The  turbulence  in  the  lake  is  modeled  by  a constant  verti- 
cal eddy  di f fusivity . This  assumption  is  almost  certainly 
invalid;  however,  lack  of  reliable  information  makes  this 
assumption  necessary  since  a satisfactory  model  for  the 
vertical  eddy  diffusivity  is  not  available. 

73.  The  displacement  of  the  lake  surface  is  assumed  to  be  small 
in  comparison  with  the  depth  of  the  lake.  This  assumption  allows  the 
complete  linearization  of  the  problem,  since  the  surface  boundary  condi- 
tions are  now  applied  at  the  undisturbed  water  surface  rather  than  the 
instantaneous  surface. 

7h.  Under  the  assumptions  listed,  the  governing  equations  are 
reduced  to  the  linear  system: 

Continuity  Equation 


Momentum  Equations 
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where  A is  the  vertical  eddy  diffusion  coefficient.  The  boundary 
conditions  are  as  follows: 


u,v,w  =0  at  z = - h(x,y) 


(67) 
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w = 0 


p = pQ(x,y)  at  z = c(x,y] 


75.  The  model  is  made  nondimens ional  by  introducing  the  dimension- 
less  variables: 
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characteristic  length  in  the  horizontal  direction 

characteristic  length  in  the  vertical  direction 

characteristic  velocity 

reference  surface  displacement 

frictional  depth 

reference  wind  shear  stress 
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W = reference  vertical  velocity 
R 

t = reference  time  value 
R 

f = Coriolis  parameters  evaluated  at  mean  latitude  of  Lake  Erie 
o 

H = mean  depth 
m 

The  nondimensional  form  of  the  governing  equations  becomes 


8u  t 3v  t 3w 
3x  3y  3z 


0 


(69) 


3u 

3t 


v 


(70) 


u 


(71) 


(72) 


p 

where  E = A /fH  is  the  vertical  Ekman  number  and  the  asterisks  have 
V v 

been  dropped  for  convenience  of  notation. 

76.  Equation  6l  can  be  integrated  over  the  vertical  to  give 


P z 


p - Po  = -z 

p = pQ(x,y)  - z (73) 

This  equation  gives  a one-to-one  correspondence  between  pressure  and 
surface  elevation;  thus,  pressure  can  be  eliminated  from  the  governing 
equations.  Normally  the  variation  of  pQ(x,y)  with  (x,y)  is  small 
enough  to  be  neglected,  Thus, 
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3p  _ 3r, 
3x  3x 


and 


ill  = & 
3y  3y 


Using  this  equation,  the  momentum  equations  in  the  horizontal  are 
written  as 


3u 

3t 


3v 

3t 


+ u = - 


- - + 

.2 

E ^ 

3x 

V3z2 

= . + 

E ^ 

9y 

V 3z2 

! 

(7b) 


(75) 

(76) 


i- 

i 

I 


and  the  vertical  momentum  equation  is  eliminated. 

II.  The  model  is  to  be  applied  to  steady-state  current  patterns; 
therefore,  the  time  dependency  in  the  governing  equations  is  eliminated. 
Thus,  basic  governing  equations  for  the  model  are  as  follows: 


^ ft+  V 


& 

dx 


E ifn  _ u = is. 

V3z2  * 


(77) 

(78) 


and 


3u  + 3v  3w 

3x  3y  3z 


0 


The  nondimensionalized  boundary  conditions  are  as  follows: 


wx 


wy 


2\ 


1 3u 
2 3z 


2X 


2 3z 


at  z = 0 


(79) 


(80) 
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w = 0 


u = v = w = 0 at  z = - h(x,y) 

where  X , a constant  which  is  proportional  to  the  ratio  of  water  depth 
to  frictional  depth,  is  given  by  X = ir  H/D. 

78.  Equations  77  and  78  can  be  combined  into  a single  equation 
by  defining  the  complex  quantities: 


t = t + i T 
w wx  wy 


K = K + i K 

3n  9x  9y 


(81) 


and  a complex  variable  representing  horizontal  velocity  T given  by 
r : U T iv 

The  complex  equation  for  F is  then  given  by: 


and 


9fr 

3z2 

ix2r 

2 

X2  9^_ 
2 9n 

conditions  are 

2 

Tw  = “2 
X 

9F 

r—  at 
9z 

z = 0 

< 

II 

O 

at 

0 

(i 

O 

II 

> 

at 

z = -h(x,y ) 

O 

II 

at 

z = -h(x,y ) 

(82) 


(83) 


This  second-order  differential  equation  involving  the  single  arbitrary 
parameter  X can  be  solved  subject  to  the  boundary  conditions  to 
yield 
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where 


u + iv 


• it  / cosh  X'z 
1 3n  ^ cosh  A 1 h 


-1 


rw  sinh  [ A 1 (h  + z )] 
2u>  cosh  ( A 1 h ) 


(81*) 


u = (1  + i)/2 

£ = surface  displacement  as  a function  of  x and  y 
A * = to  A 

h = H/H  dimensionless  depth 
m 

Z = dimensionless  coordinate 

This  equation  represents  a closed  form  solution  for  u and  v if 

3i;/3n  and  t are  known;  however,  the  surface  displacement  gradient  is 

w 

not  generally  known.  A method  must  first  be  obtained  for  evaluating 
this  quantity  prior  to  evaluating  u and  v . 

T9-  The  horizontal  mass  flux  can  be  found  by 


0 


r dz 


t a 
w 


+ h 


It 

3n 


b 


where  a and  b are  functions  only  of  the  local  depth  and  Mx  and 
are  the  x and  y components,  respectively,  of  horizontal  mass 
flux. 

8o.  The  continuity  equation  can  be  integrated  vertically  to  yield 


0 0 0 


But  w = 0 at  z 


0 and  z = -h  ; therefore. 


(85) 
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Introducing  the  mass  flux  stream  function  iji(x,y)  give; 


then 


M = ft  , M = - 
x 3y  y 3x 


(86) 


M + iM  = ? - i 
x y 3y  3x 


= t a + h b 
w 3n 


= t a + h + i b 
w ^ 3x  3y  J 


Rewriting  a = c + id  and  b = e + if  where  c , d , e , and  f are 
real  functions  of  local  depth  gives  the  following  expression  for  the 
surface  gradient: 


3£_  _ e 3£  f dip 

3x  " h(e2  + f2)  * ' h(e2  + f2)  9X 


= 

8y 


h(e2  + f2) 


dip 

9y 


h(  e£ 


f2) 


lit  + Tvx  / fc  - de  \ Twy  / c2  + df 


3x 


e2  + 


V 


2 

e + 


(88) 


The  surface  gradients  are  linear  functions  of  the  stream  function 
gradient  and  the  wind  shear.  These  equations  allow  the  surface  gradi- 
ents to  be  calculated  if  the  wind  stress  and  stream  function  are  known. 
Substituting  back  into  the  expression  for  the  complex  velocity  allows 
u and  v to  be  calculated. 

8l.  In  order  to  have  a complete  solution  for  the  problem  the 
stream  function  must  be  determined.  Equations  87  and  88  can  be  com- 
bined and  the  surface  displacement  5 eliminated  by  taking  the  3/3y 
of  the  first  equation,  - 3/3x  of  the  second  equation  and  then  adding 
the  two  equations: 
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7‘  ij;  + c 


3h 


- cr 


3h  \ 3jjj_  + 


3h 


3t 


3t 
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2 3x  + C1  3y 
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/ 3y 

3x 

+ — =*- 
3y 

3i  3t  \ 

wy  wx  j 

- n 

L 

3h 

3x  3y  / 

5 

^ wx 

3x 

wy  3y 

- c6 

A 

l vy 

3h 

3x 

3h 
wx  3y 

(89) 


where  h = Ah  and  c,  , c„  , c_  , c,  , cr  , and  are  functions 

— 1 d 5 4 p b 

of  h . This  governing  equation  for  the  stream  function  is  a linear, 
elliptic  equation  and  can  be  solved  using  the  boundary  condition  that 
the  normal  mass  flux  at  the  shore  is  zero  or  equal  to  the  river  inflow 
or  outflow, 

82.  It  should  be  noted  that  the  stream  function  formulation  is 
independent  of  the  vertical  coordinate  and  is  thus  numei ^cally  a two- 
dimensional  problem.  In  addition,  it  should  be  kept  in  mind  that  the 
stream  function  formulation  does  involve  the  vertical  eddy  diffusivity, 
which  is  basically  a function  of  the  wind  speed. 

83.  The  stream  function  equation  is  solved  numerically  by  replac- 
ing the  differential  equation  by  a finite  difference  representation. 

This  substitution  allows  the  value  of  the  stream  function  at  each  dis- 
crete point  in  the  numerical  grid  to  be  solved  in  terms  of  the  values 
at  neighboring  points  and  coefficients  that  are  functions  of  the  local 
topography  and  wind  stress  conditions.  The  resulting  set  of  simultane- 
ous equations  is  of  the  form 

[A]  M-[C] 

where  [X]  is  a square  matrix  whose  order  is  equal  to  the  number  N of 
the  interior  grid  points,  ['I']  is  the  vector  column  matrix  for  the  stream 
function,  and  [C]  is  the  vector  column  matrix  containing  the  nonhomogene- 
ous  portion  of  the  difference  equations.  The  matrix  [A]  is  sparse  and 
is  solved  using  successive  overrelaxation  by  points  (SOR)  and  by  single- 
line  overrelaxation  (SLOR). 
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8k.  Islands  in  the  lake  create  a multiply  connected  domain.  They 
can  be  included,  however,  by  requiring  that  the  value  of  ip  on  each 
island  be  specified  in  such  a manner  that  the  surface  displacement  is 
continuous  around  the  island.  This  treatment  increases  the  computa- 
tional time;  however,  islands  can  be  included  satisfactorily.  Because 
of  the  linearity  of  the  problem,  linear  interpolation  can  be  used  to 
facilitate  obtaining  solutions  for  various  wind  conditions. 

85.  Once  the  stream  function  has  been  obtained.  Equations  87  and 
88  are  used  to  evaluate  the  local  slope  of  the  lake  surface.  Once  this 
information  is  known.  Equation  8k  allows  the  local  vertical  profile  of 
horizontal  velocity  to  be  evaluated. 

86.  The  model  was  initially  applied  to  some  simple  geometries  for 
which  analytical  solutions  were  available  to  allow  the  effect  of  various 
parameters  in  the  numerical  model  to  be  considered  and  evaluated.  The 

model  was  verified  for  real-lake  conditions  by  applying  it  to  Lake 

21 

Erie.  Two  uniform  wind  conditions  were  considered  corresponding  to 
winds  of  11.8  inph  and  22.7  mph.  The  numerical  grid  step  used  in  the 
calculations  was  2 miles  in  the  open  lake  and  0.5  mile  near  the  islands 
in  the  western  basin.  The  numerical  steady-state  calculations  of  the 
wind-driven  currents  in  Lake  Erie  were  shown  to  compare  favorably  with 
current  meter  measurements  made  at  middepth  in  the  central  and  eastern 
basins.  Many  other  features  of  the  currents  observed  by  other  methods 
are  predicted  at  least  qualitatively  by  the  calculations. 

87.  The  agreement  of  calculations  and  measurements  appears  to 
indicate  that  the  shallow-lake  model  that  uses  a constant  eddy  viscosity 
is  capable  of  predicting  accurately  the  local  three-dimensional  veloci- 
ties at  middepths.  Additional  prototype  data  near  the  surface  and 
bottom  of  the  lake  in  conjunction  with  measurements  at  middepths  are 
needed  to  verify  the  model  and  the  constant  eddy  viscosity  assumption. 
The  calculated  results  show  that  for  a large  percentage  of  Lake  Erie  the 
bottom  stress  is  not  in  the  direction  of  the  horizontal  mass  flux. 

These  results  indicate  that  for  Lake  Erie  and  other  shallow  lakes  there 
may  be  large  errors  in  solutions  obtained  from  models  where  the  bottom 
stress  is  made  proportional  to  the  net  horizontal  mass  flux. 
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Three-dimensional,  time- 
dependent,  constant-density  model 
of  A.  Hag,  W.  Lick,  and  Y.  P.  Sheng 


23 

88.  The  basic  formulation  of  this  model  is  an  extension  of  the 

22 

three-dimensional  steady-state  model  of  Gedney  and  Lick.  The  method 
of  solution  of  the  model  is,  however,  completely  different.  This  model 
is  time  dependent,  and  irregular  boundaries  and  variable  bottom  topogra- 
phy can  be  handled.  It  is,  however,  restricted  to  constant-density 
flow.  The  primary  application  of  the  model  is  for  storm  surge  simu- 
lation where  surface  displacements  rather  than  velocities  are  prime 
considerations.  A free  surface  condition  is  imposed,  rather  than  the 
rigid  lid  assumption  used  in  the  steady-state  model  in  order  that  the 
details  of  the  transient  behavior  are  simulated.  As  with  the  steady- 
state  model,  the  primary  application  is  for  shallow  lakes  with  param- 
eters in  the  program  closely  corresponding  to  those  for  Lake  Erie. 

89.  Equations  63-66  and  the  boundary  conditions  given  by  Equa- 
tions 67  '..and  68  form  the  basic  system  of  equations  fot*  the  model.  These 
equations  are'  put  in  a nondimensional  form  with  some  slight  changes  from 
the  nondimensionalization  process  used  by  Gedney  and  Lick.  The  non- 
dimensional  variables  are  introduced  as  follows: 


x*  - — v*  = X-  z#  = — 

L ’ y L ’ ^ H 


# - H__  # - X_  * _ Lw 

U = UR  ’ " = UR  ’ " = HUR 


,*  = — £ — r*  = t*  = f+ 

pflLL  ’ ; fU  L » L 


(90) 
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90.  The  nondimensionalized  forms  of  the  governing  equations 


3u  + 3v_  + 3w  _ 
3x  3y  3z 


8iL_v  = _i£+E; 

3t  3x  V „ 2 

3z 


_ i£+  . 9_v 
3y  N ,2 


in  , kh 

3z  ' fUDL 

R 


where  the  asterisks  have  been  dropped  for  convenience. 

91.  The  dimensionless  boundary  conditions  are: 


u,v,w  = 0 at  z = -I  — 


3z  x 


at  z = 0 


3v 

3z  “ Ty 


at  z = 0 


2 2 

f L 3t 

W = — TT  at  Z = 0 

gH  3t 


92.  The  hydrostatic  Equation  9^  is  integrated  over  the  lake  depth 


to  obtain 


P = PQ(x,y)  ♦ C - z 

R 
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Assuming  that  the  atmospheric  pressure  is  uniform  results  in 


»£  = K 

3x  3x 


3£  _ 3^ 
3y  3y 


Substituting  into  Equations  92  and  93  results  in 


= _ i£+  E,  Ui 

3X  ^ 3z2 


+ u = _ 

3t  3y  .2 


93.  Equation  91  is  further  modified  by  defining 


U = I u dz 


V = | v dz 


The  continuity  equation  is  then  integrated  over  the  depth  of  the  lake 
to  obtain 
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o 


o 


o 


For  nondimens ional  variables  the  continuity  equation  becomes 


f?L2  + 9U  3V  _ 
gH  3t  3x  dy 


or 


+ /in  + il) 

3t  f2L2  ^3x  3y/ 


(102) 


9*+.  Equations  98,  99,  and  102  must  be  solved  to  obtain  the  veloc- 
ity field  and  surface  elevations,  £ , in  the  lake  using  this  model. 
These  equations  are  subjected  to  the  boundary  conditions  given  by 
Equation  95.  The  model  as  formulated  does  not  provide  for  a mass  flux 
inflow  or  outflow  through  the  boundaries  of  the  lake. 

95-  To  insure  no  loss  of  accuracy  in  the  numerical  computation  in 
the  shallow  regions  of  the  lake,  a more  convenient  coordinate  system  is 
introduced.  The  desired  transformation  maps  the  bottom  of  the  lake  on- 
to a constant  a surface  where  o is  a transformed  vertical  coordi- 
nate. A uniform  o grid  is  then  used  in  the  (x,y,o)  system  that 
corresponds  to  a nonuni  form  vertical  grid  in  the  (x,y,z)  system.  The 
new  a system  is  obtained  by  replacing  the  vertical  coordinate  in  the 
(x,y,z)  system  by 
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/ 


0 = 1 + 


h(x,y ) 


where  h(x,y)  is  the  dimensionless  local  depth. 

96.  The  governing  equations  can  then  be  written  as 


3u 

3t 


F p 

ac  , v a u 

ax  , 2 ? 

h 3o 


where 


& + u = _ & + h A 

3t  3y  h2  ao2 


2 2 
L f 


and  the  horizontal  mass  fluxes  U and  V are  defined  by 


U = h(x,y ) J u da 

o 

1 

V = h( x,y ) ^ v da 
o 


The  boundary  conditions  are 
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(103) 


(1010 


(105) 


(106) 


(107) 


X 


u = v = 0 at  a=0 


3u  , 

r — = hT 

3o  x 


at  a = l 


(108) 


3v  , 
t—  = hx 
3o  y 


at  o = l 


and  V = 0 at  the  shoreline  where  V is  velocity  normal  to  shore- 
n n 

line. 

9f.  The  finite  difference  grid  is  set  up  using  constant  grid 
spacings  Ax  , Ay  , and  A o and  generally  Ax  equals  Ay  . The 
variables  are  defined  as  shown  in  Figure  2b.  The  staggering  of  the 
variables  simplifies  the  representation  of  the  required  derivatives. 
Centered  differences  in  space  and  forward  differences  in  time  are  used 
in  the  numerical  calculations.  Values  for  variables  are  sometimes 
needed  at  points  at  which  they  are  not  defined.  These  values  are 
obtained  by  linearly  interpolating  from  the  values  immediately  adjacent 
to  the  point  under  consideration. 

98.  Equations  103-105  are  written  in  a finite  difference  form 
using  forward  time  and  centered  space  differences  and  for  stability 
consideration  a modified  DuFort-Frankel  scheme  for  the  vertical  dif- 
fusion terms.  The  formulation  is  explicit.  In  addition,  the  mass 
fluxes  U and  V are  evaluated  by  numerically  integrating  Equa- 
tions 106  and  107  using  Simpson's  rule.  The  finite  difference  equations 

at  time  t , can  then  be  written  as 
n+1 


n+1 

Ui,j+l/2,k 


1 

x + a 


U i , j+l/2,k  + At 


r n _ fd^\n 

(_V<i,j  ,+l/2,k>  " \dxj  <i  , J+l/2> 


— fun 

At  ^ i.J+1/2, 


n-1  n 

k+1  " i , j+1/2  ,k  Ui,j+l/2 


for  1 < k < KM( surface) 


(109) 
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LEGEND 

X POINTS  DEFINING  U,  u (DEFINED  AT  i,  Jt'/2,K) 

D POINTS  DEFINING  V,  v (DEFINED  AT  i + l/2,  j,  k) 

O POINTS  DEFINING  5 (DEFINED  AT  i ♦ '/2 , J + K) 

V POINTS  DEFINING  u 

A POINTS  DEFINING  V 

Figure  2h . Spatial  grid  used  in  the  numerical  computations 

8U 


n+1 

Ui.,,j+1/2,KM 


1 + a 


,J+1/2,KM 


+ At 


n 

V<i ,J+l/2,KM> 


<i,J+l/2> 


+ 


— ( ?un 

At  \ i , j+l/2,KM-l 


n-1 

Ui,j+1/2,KM  + 


(no) 


V 


n+1 

i+l/2  ,j 


,k 


n 

u<i+l/2, j ,k> 


.J> 


a [ n 
At  \^Vi 


i+l/2, J,k+1 


V / + 

i+l/2, J,k 


i+l/2, J,k-1 


(111) 


for  1 < k < KM  (surface) 
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where 

a = EyAt/lh'A  o2) 

KM  corresponds  to  o = 1 (surface) 


Y = rr/(pURHf) 
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points  at  which  variables  must  be  interpreted  since  they  are 
not  evaluated  in  the  staggered  grid 


99-  Thus,  the  numerical  procedure  can  be  summarized  as  follows: 

a.  Assume  that  the  values  at  time  t are  known  (initially 
at  t,  u=v=?=0). 

b.  Use  Equations  108  and  109  to  calculate  un+^  . 

c_.  Use  Equations  111  and  112  to  calculate  vn+1  . 

d_.  Evaluate  Un  ^ and  by  numerically  integrating 

Equations  106  and  107. 

£.  Set  U and  V equal  to  zero  on  appropriate  boundary 
segments . 

_f.  Evaluate  c;n+1  from  Equation  113. 

The  calculation  is  completed  at  time  t and  the  computation  for  the 
next  time  step  can  now  be  started  with  step  a. 

100.  The  numerical  model  consists  of  a system  of  five  differential 
equations  together  with  two  numerical  integrations.  The  stability  of 
the  system  is  found  to  be  limited  only  by  the  well  known  limitation  on 
surface  gravity  waves. 

101.  The  model  was  initially  verified  by  applying  it  to  some  simple 

lake  geometries:  an  infinitely  long  lake,  a square  basin  with  constant 

depth,  and  a constant  bottom  slope  basin.  The  numerical  solutions  were 
compared  with  analytical  solutions  to  determine  the  effects  of  changes 
in  the  eddy  diffusivity,  lake  geometry,  and  bottom  topography. 

2h 

102.  The  Wilson  formulation  for  the  wind  stress  is  used  in  the 
model : 


T 

w 


p C,  W W 
ad1  a 1 a 


CllU ) 


where  p is  the  density  of  the  air,  C is  the  drag  coefficient  (C  = 

3.  Q Q 

0.00273  for  winds  > 6 m/sec  and  = 0.00166  for  winds  < 6 m/sec),  and 

W is  the  wind  velocity  at  10  m above  the  water  surface.  For  a 
a 

spatially  varying  and  time-varying  wind  stress,  the  stress  at  any  point 
in  the  lake  is  obtained  by  interpolating  from  the  known  value  at  a few 
velocity  measuring  stations  around  the  lake.  The  method  of 

86 

* 


interpolation  is  similar  to  that  used  by  Platzman. 


lU 

103.  The  numerical  model  has  been  applied  to  Lake  Erie  for  two 
22 

cases.  The  first  case  is  for  a uniform  wind  blowing  along  the  long 
axis  of  the  lake.  The  actual  verification  of  the  model  is  obtained  by 
calculating  the  flow  in  the  lake  during  storm  Agnes.  The  calculated  and 
measured  values  for  the  surface  displacement  at  Cleveland  and  at  Port 
Stanley  are  found  to  be  in  good  agreement.  Also,  the  spatial  variation 
of  the  surface  displacement  is  consistent  with  the  limited  amount  of 
data  that  are  available. 

lOU.  In  applying  the  model  to  Lake  Erie,  it  should  be  noted  that 
much  of  the  western  basin  and  regions  very  close  to  shore  become  very 
shallow.  In  such  regions,  the  assumption  that  the  surface  displacement 
is  small  compared  with  the  depth  may  be  inadequate,  and  the  linearized 
theory  is  likely  to  be  inaccurate.  The  minimum  depth  considered  in  the 
numerical  model  was  1.8  m.  In  addition,  the  islands  separating  the 
western  and  central  basins  of  Lake  Erie  are  neglected  in  the  present 
model. 

105.  A uniform  horizontal  grid  with  x = y = 6.U  km  is  used  in 
the  entire  lake  for  the  Lake  Erie  calculations.  In  the  vertical  direc- 
tion, the  region  from  the  top  to  the  bottom  (0  = 0 to  0 = l)  is 
divided  into  five  equal  intervals  of  0 . In  actual  physical  variables, 
the  lengths  of  these  intervals  are  smaller  in  shallow  regions  and  larger 
in  the  deep  regions  of  the  lake.  The  time  step  used  in  the  calculations 
for  Lake  Erie  was  50  sec,  and  the  results  showed  no  evidence  of 
instability. 

106.  Calculations  were  made  for  eddy  diffusivity  values  of 

2 2 2 
25  cm  /sec  and  40  cm  /sec.  For  A^  = 25  cm  /sec  the  numerical  results 

for  surface  displacement  as  a function  of  time  for  storm  Agnes  compare 
favorably  with  the  measured  values  of  the  surface  displacement.  The 
surface  displacement  values  for  A^  = hO  cm  /sec  are  indicated  as  being 
practically  indistinguishable  from  those  for  A^  = 25  cmc/sec  indi- 
cating that  the  magnitude  and  direction  of  the  wind  velocity  dominate 
the  flow.  This  situation  would  not  exist  after  the  storm,  as  the  sur- 
face oscillates  and  decays  to  an  equilibrium  position.  The  accuracy 
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of  velocities  predicted  by  the  model  cannot  be  properly  evaluated 
because  of  a lack  of  appropriate  prototype  data. 

107.  The  authors  reference  the  numerical  model  study  of  storm 

20 

Agnes  on  Lake  Erie  by  Paskausky,  whose  model  is  vertically  integrated 
and  contains  a rigid  lid  assumption.  Paskausky' s values  of  setup  at 
Cleveland  and  setdown  at  Port  Stanley  are  indicated  as  being  con- 
sistently lower  than  the  measured  values.  The  authors  attribute  this 
difference  to  the  use  of  the  rigid  lid  assumption  combined  with  a 
greater  minimum  allowable  depth  (2b  ft)  in  Paskausky 's  model.  To  verify 
their  conclusions,  the  authors  made  several  calculations  with  a verti- 
cally integrated  model  using  a minimum  allowable  depth  of  6 ft.  This 
model  predicted  much  better  agreement  with  measured  values  than  did  that 
of  Paskausky’ s. 

Three-dimensional  variable- 
density  model  by  Jan  J.  Leendertse, 

Richard  C.  Alexander,  and  Shiao-Kung  Liu 

25 

108.  This  model  lays  the  foundation  for  the  development  of  a 
layered  model  to  be  used  for  numerical  simulation  of  the  fluid  flow  in 
water  bodies  with  irregular  boundaries  and  nonisotropic  density  distri- 
bution. Work  on  the  model  is  not  complete,  and  only  very  limited  appli- 
cations and  attempts  at  model  verification  are  reported.  The  model  is 
intended  for  application  to  estuaries  and  coastal  seas.  Some  of  the 
assumptions  and  formulation  procedures  are  not  immediately  adaptable 

to  freshwater  lakes,  where  density  variations  are  due  to  temperature 
variations  rather  than  to  salinity  variations.  The  general  formulation 
procedures  can  be  extended,  however,  to  be  applicable  to  freshwater 
lakes . 

109*  The  conservation  form  of  the  governing  equations  is  used  to 
allow  the  preservation  of  certain  integral  conservative  relations  of 
the  continuum  equation  in  the  finite  difference  formulation  of  the 
problem.  This  form  insures  that  mass,  momentum,  or  other  variables  are 
neither  created  nor  destroyed  as  a result  of  the  computational  scheme. 
The  fluid  density  is  assumed  to  vary  with  salinity.  The  basic  governing 
equations  for  the  formulation  are  those  given  by  Equations  59-62. 
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110.  In  addition,  because  the  density  is  considered  ar,  a function 
of  salinity,  S,  a constitutive  equation  and  salinity  transpnr  . and 
diffusion  equation  are  required. 

Salinity  transport  and  diffusion 


Constitutive  equation 


p = p + p ' (S) 


(116) 


where 


D ,D  ,D 
x’  y’  z 


P 

p' 


= salinity 

= eddy  diffusion  coefficients  of  salinity  in 
and  z directions 

= constant  reference  density 

= departure  from  p depending  on  salinity 


x 


y , 


111.  For  a more  general  constitutive  equation,  a density  depend- 
ence on  temperature  as  well  as  on  salinity  would  be  required.  This 
more  general  formulation  would  require  an  additional  equation  governing 
the  energy  balance  for  the  system.  Diffusion  equations  for  the  various 
other  substances  of  interest  would  also  be  required.  The  actual  numeri- 
cal simulation  of  the  three-dimensional  flow  problem  in  this  model  is 
considerably  more  complex  than  for  the  previous  models  that  have  been 
considered.  Fewer  simplifications  of  the  basic  equations  are  made,  and 
the  resulting  model  is  highly  nonlinear  with  the  accompanying  diffi- 
culties of  formulating  the  nonlinear  terms. 

112.  The  basic  approach  in  the  layered  model  is  to  visualize  the 
fluid  motion  in  horizontal  slices.  Within  each  horizontal  slice  the 
equations  are  integrated  over  the  height  of  the  layer  to  obtain  at  any 
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horizontal  position  depth-averaged  fluid  properties  for  the  slice.  For 
any  slice,  this  approach  is  essentially  identical  with  the  two- 
dimensional  depth-averaged  model  formulation.  The  layers  are  not  inde- 
pendent, however,  since  there  is  an  exchange  of  mass,  momentum,  and 
salinity  between  layers.  The  air-water  interface  is  the  boundary  for 
the  system's  upper  layer,  and  the  bottom  defines  the  boundary  condition 
for  the  bottom  layer.  At  these  boundaries  the  mass  and  salinity  fluxes 
are  zero.  Only  momentum  is  transferred,  at  the  surface  as  a driving 
force  by  wind  and  at  the  bottom  as  a dissipative  effect  due  to  friction. 

113.  Figure  25  shows  a sample  layout  of  the  vertical  grid.  The 
surfaces  z = z = constant  are  levels  separating  the  various  layers 

K 

of  thickness  h . The  thickness  h.  will  generally  vary  in  space  and 

K.  -L 

time  due  to  wave  action.  The  intermediate  thicknesses  are  constant  over 
the  region  while  the  bottom  layer  will  vary  with  x and  y according 
to  the  prescribed  bottom  topography.  The  number  of  layers  will  vary 
from  one  position  in  the  horizontal  plane  to  another  depending  on  the 
depth. 

llU.  The  governing  equations  are  integrated  over  the  thickness 
of  the  layer  in  the  same  manner  as  was  worked  out  in  detail  previously 
for  the  two-dimensional  depth-averaged  models.  The  momentum  equations 
for  the  k layer  become 


k-1/2  k-l/2  k-1/2  k-1/2 


k+l/2  k+1/2  k+l/2  k+l/2 
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The  equations  can  be  written  in  the  form: 
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(117) 
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where  is  the  depth-averaged  density  for  the  K layer. 

115-  Some  approximations  are  necessary  at  this  point.  The  second 
and  third  terms  represent  integrals  of  products  of  velocities.  These 
terms  are  replaced  by  products  of  integrals. 


_3 

3x 


<uu> 


k-1/2 


k-1/2 


k-1/2 


uu  dz 


_3 

3x 


k+1/2 


k+1/2 


u dz 


u dz 


k+1/2 


Similarly 


and 


3(uv) 3 (l 

3y  ~ 3y  \h 


<u> 


3 ( uv) 3/l 

3x  ~ 3x  \h 


<u> 


3(w) 

3y 


_3 

3y 


& 


<v> 


<v> 


(118) 


The  density  is  considered  to  vary  only  with  horizontal  position  wi t.h i n 
each  layer.  The  differentiation  and  integration  processes  are  taken 
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as  interchangeable  in  the  momentum  diffusion  terms  and  are  approximated 
by: 
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116.  In  integrating  the  hydrostatic  pressure  variation  across  the 

layer 


k-1/2  k-l/2 


k+1/2  k+1/2 


approximations  must  be  introduced  since  p and  p are  not  defined  at 
the  interface  between  elements  but  at  the  center  of  the  elements. 
Therefore 


k-1/2 


k+1/2 


iti 

9z 


dz 


can  be  approximated  by: 


p - p 
k k-1 


where  ?k  is  the  layer-averaged  pressure  of  the  kth  layer  and 

k-1/2 

pg  dz 

k+1/2 

is  approximated  by 
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The  same  principle  is  extended  to  represent  all  the  interface  values  of 
dependent  variables  as  the  two  point  averages  of  the  adjacent  layer 
values  for  the  same  variable;  for  example, 

Uk+l/2  = 2 (Vk  + Vk+l)  (120) 

where  is  the  layer-averaged  U-component  of  velocity. 

117.  From  the  hydrostatic  pressure  variation  an  expression  for  p 
can  be  obtained: 


Pg 
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or  simply 


pg  dz  + f ( x 
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A 


if  the  variation  in  atmospheric  pressure  over  the  water  body  is  neg- 
lected. A change  in  pressure  at  a particular  location. can  therefore 
occur  as  a result  of  an  effective  change  in  depth  (change  in  t ) or  as 
a result  of  a variation  in  density. 

118.  The  horizontal  pressure  gradient  term  for  the  upper  layer 
can  then  be  approximated  by: 
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(121) 


The  pressure  gradient  for  the  other  layers  is  given  by 
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(122) 


where  the  total  number  of  layers  B is  a function  of  local  depth. 

119.  When  the  continuity  equation  is  integrated  over  layer  k , 
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Interchanging  the  differential  and  integral  operators 
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for  a particular  layer  gives 
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(123) 


At  the  bottom  w = 0 ; therefore 
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At  the  water  surface  w = 9^/3t  and  therefore 


1=1 


9 ( h1v1 ; 


= 0 


(125) 


reflects  the  free-surface  conditions. 

120.  Horizontal  momentum  diffusion  terms  are  expressed  as  gradi- 
ents of  Reynolds  stresses  evaluated  by  means  of  eddy  viscosity  coef- 
ficients 
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. 3v 
T = A — 
yx  yx  3x 


These  terms  are  further  approximated  hy  assuming  only  a single  eddy 
viscosity  coefficient  characteristic  of  the  u and  v components  of 
velocity,  i.e. , 
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121.  The  final  governing  equations  for  the  layer  are 
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where 
stress . 


x and  t are  the  x-  and  y-  components  of  inter facial  shear 
xz  yz 

At  the  bottom 
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and  at  the  surface 
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122.  The  boundary  stresses  are  handled  in  a conventional  manner. 

At  the  bottom  the  stress  is  related  to  the  Chezy  coefficient  C . 

/ 2 A 2 
Pg  guy  u + v 

TBx=  ^ 

(133) 

/ 2 A 2 
pb  gvyu  + v 

TBy  ^ 


At  the  surface  the  shear  stress  is  formulated  in  terms  of  the  quadratic 

law. 
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2 

T = pC  U 
d 


(131* ) 


where  u = mean  wind  speed 

123.  The  equation  of  state  is  taken  as 


P 


Po 


VVo 


(135) 


While  X'  can  be  expressed  as  a function  of  both  temperature  and  sa- 
linity, p is  considered  to  vary  only  with  salinity  within  this  par- 
ticular problem  formulation. 

12l*.  A space-staggered  grid  was  selected  for  the  finite  difference 
representation  of  Equations  128-132.  The  grid  structure  is  shown  in 
Figures  26  and  27. 

125.  Because  of  the  nonlinear  nature  of  the  problem,  the  finite 
difference  formulation  is  complex  and  cumbersome.  Sum  and  difference 
operations  are  introduced  as  an  aid  in  formulating  the  finite  difference 
scheme.  For  an  arbitrary  variable 

F H F (iAx,  JAy,  kAz,  nAt ) 

The  following  sum  and  difference  operators  are  adapted  for  representing 
the  function  at  (i,  j,  k,  n)  : 


FX  = ~ ^F[ ( i + |-)Ax,  JAy,  kAz,  nAt] 

+ F[ ( i - Ax,  JAy,  kAz,  nAt]j  (136) 

SXF  = + |-)Ax,  JAy,  kAz,  nAt] 

- F[(i  - |-)Ax,  JAy,  kAz,  nAt]|  (137) 
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Figure  26.  The  location  of  u (-),  v ( 1 ) , and  other 
parameters  (+)  in  the  space-staggered  grid 
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Figure  27 . Relative  position  of  the  variables  in  the  model 
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Similar  operations  are  defined  for  y , z and  t . The  special 
notation  used  for  shifting  time  levels  is 

F+  = F [i  Ax  , j Ay  , k Az  , (n+l)  At] 

F_  = F fi  Ax  , j Ay  , k Az  , (n-l)  At]  (138) 

l:’(..  The  finite  difference  approximations  of  Equations  128-132 
used  for  integration  become 

= - ? + «y(hyv)]  (139) 

at  i , j , n 
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6 (hS)1  = - 6 (hXuSX)  - 6 (hyvSy)  - hZ6  (vSZ) 
t x y z 


+ hZ6  ( k6  S) 
z z 

at  i,  j,  k,  n 


(1U2) 


These  finite  difference  equations  in  the  compact  notation,  which  was 
used,  have  the  same  appearance  as  the  differential  equations. 

127.  ' The  density  is  computed  with  the  equation  of  state  as  follows: 


P 


, A + ix  p 
o o 


at  i,  j,  k,  n + 1 


(11*3) 


128.  The  following  finite  difference  equations  are  used  to  compute 
derived  variables: 


6 w = - 6(hXu)  - 6(hyv) 
z 

at  i,  j , k,  n + 1 ( lU  U ) 

6 p = g(pX6  5)  + 1/2  (hX6  p) 

X X X 

at  i + 1/2,  j,  1,  n + 1 ( 1 U 5 ) 

6 p = g(py6  c)  + 1/2  (hy6  p) 

J J J 

at  i,  j + 1/2,  1,  n + 1 (lU6) 

6 (6  p)  - g hZS  pZ 
z x x 

at  i + 1/2,  j,  k + 1/2,  n + 1 (1U7) 


6z(6yp)  = g hZ6ypZ 

at  i,  j + 1/2,  k + 1/2,  n + 1 ( lU8 ) 
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Boundary  stress  and  inter- 
facial  stress  term  formulations 


129*  At  the  surface,  the  stress  term  can  be  computed  from 


( txz\  = C*  — w2  sin  4>  at  i + 1/2,  j,  1,  n (ll*9) 

VpXZ  " / k-1/2  pX  a 


( ~~T  xYZ ) = c*  vn  cos  <t>  at  i.  J + 1/2,  1,  n (150) 

\pyZ  " /k-1/2  py 


where 

C*  = resistance  coefficient 

p^  = atmospheric  density 

w = wind  speed  at  10-m  level 
a 

<{)  = angle  between  wind  direction  and  y-axis 
Similarly,  the  stress  term  at  the  bottom  (in  layer  B)  can  be  computed 
from  Equation  133. 
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at  i,  j + 1/2,  B,  n 


130.  Since  estuary  motions  are  mostly  turbulent,  the  effects  of 
eddy  viscosity  play  a significant  part  in  controlling  stresses  between 
adjacent  layers  with  different  velocities.  Mass  and  momentum  transfers 
are  the  result,  and  these  transfers  tend  to  equalize  relative  velocity 
differences  in  both  the  lateral  and  the  vertical  directions.  If  the 
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quadratic  relationship  between  the  interlayer  stresses  and  the  velocity 
differences  is  assumed  applicable,  the  following  expressions  for  the 
interlayer  stresses  can  be  used. 
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(153) 


(15*+) 


where  v is  the  coefficient  of  interfacial  friction. 

131.  The  equations  used  in  the  numerical  integrations  are  second- 
order  approximations  in  time  and  space,  except  that  the  evaluation  of 
the  momentum  diffusion  by  use  of  Equations  153  and  15*+  in  Equations  l*+0 
and  1*+1  is  not  central  in  time,  but  at  a lower  time  level,  since  other- 
wise the  computation  would  become  unstable.  For  the  same  reason  the 
dispersion  terms  in  Equation  l*+2  are  taken  at  the  lower  time  level. 

132.  In  describing  the  progress  of  the  computation,  it  is  assumed 
that  the  computation  has  progressed  to  a time  level  n . At  this  time 
the  velocity  components  u and  v , the  salinity  S , and  the  pressure 
p are  all  determined  for  each  point  of  the  grid  and  available  for 
computation.  These  variables  are  also  available  at  the  previous  time 
step  n - 1 . The  water  levels  are  available  for  computation  at  time 
level  n and  n - 1 . 

133.  At  the  boundaries  of  the  water  body  to  be  computed,  all 
diffusion  coefficients  are  zero,  as  are  the  velocities  perpendicular 

to  the  boundary.  In  this  manner,  no  mass  fluxes  or  diffusive  transports 
of  salt  will  result.  The  boundaries  of  the  horizontal  layers  also 
satisfy  a no-slip  boundary  condition. 

13*+.  With  this  information,  the  prediction  equations  139-l*+2  are 
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used  to  advance  the  computational  field  of  £ , u , v , and  S . 

For  the  evaluation  of  the  boundary  and  interfacial  stress  terms  in  these 
equations.  Equations  1U9— 15^»  are  used.  These  are  all  explicit  oper- 
ations for  each  of  the  variables  resulting  in  computed  values  of  r,  , 
u , v , and  S at  time  level  n + 1 . 

135.  Subsequently,  the  vertical  velocities,  the  densities,  and  the 
pressures  required  for  the  next  step  of  the  integration  are  computed. 

For  the  vertical  velocity.  Equation  lU 3 is  used,  starting  at  the  bottom 
layer  and  proceeding  upward  with  the  computation.  In  this  manner  the 
vertical  velocity  at  the  lower  interface  is  known,  and  the  other  verti- 
cal velocities  can  be  determined  explicitly. 

136.  The  density  for  each  grid  point  can  be  computed  by  using 
Equation  135  from  the  salinity  data  now  available  at  time  n + 1 . As 
water  levels  and  densities  are  known  at  time  level  n + 1 , the  pressure 
gradients  in  the  x and  y directions  can  now  be  computed  by  using 
Equations  lH5-ll+8.  In  this  case,  the  gradients  at  the  surface  layer, 
and  subsequently  the  other  layers,  proceeding  downward,  are  computed 

by  using  Equations  1^5  and  lU8.  In  this  manner  the  horizontal  pressure 
gradient  in  the  upper  layer  of  the  two  layers  involved  in  Equations  1^7 
and  1^8  is  known,  and  the  horizontal  pressure  gradient  in  the  lower 
layer  can  be  determined  explicitly. 

137*  The  computational  cycle  is  now  ready  to  be  reported,  and  c,  , 
u , v , and  S can  be  computed  for  time  level  n + 2 , as  all  vari- 
ables are  now  determined  for  time  level  n + 1 . To  start  a compu- 
tation, information  at  two  time  levels  is  required. 

138.  The  model  was  initially  applied  to  several  simple  cases  to 
test  the  computational  procedure.  These  test  cases  were  characterized 
by  an  increasing  complexity  in  both  boundary  geometry  and  hydrodynamic 
behavior  of  the  process.  These  test  cases  included 

a.  A seiche  oscillation  in  a rectangular  basin  (density 
constant ) . 

b.  Wind-driven  circulation  in  a rectangular  basin  with  wind 
direction  along  the  axis  of  the  rectangle  (density 
constant ) . 
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c_.  Wind-driven  circulation  is  a rectangular  basin  with  a 
diagonal  wind  stress  (density  constant). 

d_.  Wind-driven  circulation  in  a rectangular  basin  with  a 
density  gradient. 

e_.  Wind-driven  circulation  in  Lake  Michigan  with  a constant 
density. 

The  model  results  for  the  test  case  appear  reasonable  and  verify  the 
general  operating  principles  of  the  model.  No  attempt  has  been  made 
however  to  simulate  real  events  using  observed  field  data.  It  should 
also  be  noted  that  no  model  results  are  indicated  for  the  general  prob- 
lem of  variable  geometry,  variable  topography,  and  variable  density. 

Work  on  this  model  is  still  progressing.  To  make  the  model  applicable 
to  engineering  investigations,  the  following  development  work  is  recom- 
mended : 

a.  Formulate  and  develop  the  computational  code  for  bounda- 
ries, such  as  those  that  occur  at  the  bottom,  sides,  and 
the  seaward  ends  of  the  model,  based  upon  realistic 
formulations  of  the  local  physical  processes  and 
conditions . 

b.  Investigate  the  transfer  of  energy  from  one  frequency 
range  to  another  in  the  model  by  analyzing  properties  of 
the  finite  difference  scheme  and  by  experimenting  with 
different  formulations  of  the  advective  terms  in  the 
equations  of  motion. 

c_.  Develop  methods  for  graphically  representing  results. 

d_.  Incorporate  the  transport  of  heat  in  the  model,  together 
with  the  interaction  of  temperature  distribution  with 
the  fluid  flow. 

139-  This  model  presently  provides  for  a variation  of  density  only 
with  salinity  and  thus  is  not  directly  applicable  to  Lake  Erie  as  a 
variable  density  model.  It  could  be  applied  to  Lake  Erie  as  a three- 
dimensional  constant  density  model.  The  other  alternative,  as  far  as 
application  to  Lake  Erie,  is  to  introduce  the  energy  equation  and  con- 
vert the  model  to  treating  density  as  a function  of  temperature  if  a 
variable  density  lake  model  is  desired.  The  author  indicates  that  this 
modification  is  presently  being  investigated. 
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Model  of  J.  Paul, 

W.  Lick, and  R.  T.  Gedney 


lUO.  A t ime-dependent , three-dimensional,  viscous,  variable- 
density,  heat-conducting  model  is  presently  being  developed  by  J.  Paul 
and  W,  Lick  of  Case  Western  Reserve  University  and  R.  T.  Gedney  at 
National  Aeronautics  and  Space  Administration  (NASA)  Lewis  Research 
Center.  This  model  is  not  yet  completely  operational;  however,  it  is 
an  extension  of  an  operational  model  for  thermal  plumes  and  river  dis- 
charges  developed  by  Paul  and  Lick/'  The  thermal  plume  and  river  dis- 
charge model  was  developed  for  studying  the  near  field  surrounding  the 
point  of  discharge  of  a river  or  power  plant  discharging  into  a body  of 
water.  Application  of  the  thermal  plume  model  required  establishing 
computational  limits  in  the  body  of  water  to  which  the  plume  is  being 
discharged  and  appropriate  boundary  conditions  on  the  computational 
boundaries.  In  the  new  model,  the  plume  model  will  be  coupled  with  a 
variable  density  model  for  a freshwater  lake  to  allow  the  entire  lake 
and  ihe  plume  to  interact.  Changing  lake  conditions  will  then  be  re- 
flected in  an  appropriate  change  in  the  plume. 

lUl.  The  formulation  for  the  lake  model  is  not  available  at 
present;  it  must,  however,  be  consistent  with  the  thermal  plume  model 
where  the  two  solutions  are  meshed  together.  The  actual  model  formu- 
lation is  more  in /olved  than  any  model  discussed  previously  in  this 
report  due  to  the  introduction  of  an  energy  balance  equation. 

lU2.  While  the  complete  formulation  of  the  model  is  not  available 
at  present,  the  plume  model  will  be  discussed  as  illustrating  some  of 
the  methods  and  techniques  required  in  the  model.  The  governing  equa- 
tions for  the  plume  model  are  derived  from  the  time-dependent,  three- 
dimensional  equations  of  motion  for  a viscous,  compressible,  heat- 
conducting  fluid  with  a linear  stress-strain  relationship. 

1J*3.  The  time-dependent,  three-dimensional  equations  of  motion 
for  a viscous,  compressible,  heat-conduting  fluid  with  a linear  stress- 
strain  relationship  are  as  follows: 
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density 

absolute  viscosity 

specific  heat  at  constant  pressure 

absolute  temperature 
thermal  conductivity 
heat  dissipation 
heat  source  and/or  sink 
equation  of  state 


Note  that 
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where 


The  density  is  taken  to  be  dependent  only  upon  the  temperature. 

The  geometry  of  the  plume  is  shown  in  Figure  28  where  b 

o 

and  hQ  are  width  and  depth,  respectively,  of  river  discharging  the 

plume.  In  obtaining  the  governing  equations  for  the  plume  model,  the 

following  assumptions  are  made: 

a.  The  pressure  is  assumed  to  vary  hydrostatically,  and 
therefore 


b.  The  rigid  lid  approximation  is  made,  i.e.,  w(z=0)  = 0 . 

c_.  The  Boussinesq  approximation  is  made.  This  approximation 
assumes  that  the  density  variations  are  small  and  can  be 
neglected  except  in  the  gravity  term. 

d_.  Heat  sources  and/or  sinks  in  the  fluid  are  neglected. 

e_.  Eddy  coefficients  are  used  to  account  for  the  turbulent 
and  molecular  diffusion  effects  in  both  the  momentum  and 
energy  equations.  The  horizontal  coefficient  is  assumed 
to  be  constant,  but  the  vertical  coefficient  is  assumed 
to  be  dependent  on  the  local  vertical  temperature  gradient. 

_f.  The  variations  in  the  bottom  topography  are  assumed  to 
be  gradual. 


All  of  these  assumptions  except  d_  and  _f  along  with  their  consequences 
have  been  discussed  previously  in  connection  with  other  models.  All 
heat  inputs  and  outputs  to  the  model  are  assumed  to  occur  at  the 
boundaries  of  the  model.  As  a consequence,  heat  transfer  by  radiation 
to  the  water  is  treated  as  a surface  heat  flux. 
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, 27 

lU^.  Paul  and  Lick  presented  in  a previous  study  the  following 

detailed  mathematical  description  of  their  plume  model. 

The  present  numerical  model  allows  variations  in 

the  depth  of  the  basin  which  the  outfall  discharges 

into.  A standard  numerical  procedure  to  fit  the 

variable  depth  into  a model  is  to  vary  the  number  of 

vertical  points  in  the  computational  mesh  according 
r po  I 

to  the  local  depth.  c A seemingly  more  complicated 
procedure,  although  a lot  simpler  in  many  aspects, 
is  to  stretch  the  vertical  coordinate  with  respect  to 
the  local  depth.  The  equations  are  transformed 
according  to 

x *-*■  x, 

y *-*  y, 

a *-*■  z/h(x,y). 

The  equations  to  be  solved  are  more  complicated  look- 
ing because  of  the  appearance  of  the  depth  in  the 
equation,  but  they  are  solved  for  a basin  of  constant 
depth  in  the  transformed  system.  This  greatly  reduces 
the  programming  complexities  of  the  model  and  makes 
the  inclusion  of  depth  variations  simpler. 

The  assumption  of  gradual  variations  in  the 
depth  allows  a reduced  form  for  the  transformed 
diffusion  terms  to  be  used.  The  transformation  used 
is  not  conformal  and  so  the  transformed  diffusion 
terms  involve  cross-derivatives  of  the  spatial 
coordinates.  The  terms  containing  derivatives  of 
the  depth  are  neglected  with  respect  to  those  terms 
containing  only  the  depth.  This  approximation  is 

used  in  meteorological  problems  when  topographic 

[29  30 ] 

variations  are  included.  ’ 

The  resulting  system  of  transformed  equations, 
in  nondimens ional  form,  are  the  following: 
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p = p(T  ) and 

O Cj 

uq  = reference  velocity, 

b = horizontal  reference  length, 
o 

h = vertical  reference  length, 
o 

Ajj  = horizontal  eddy  viscosity, 

Ay  = vertical  eddy  viscosity, 

B„  = horizontal  eddy  diffusivity. 
By  = vertical  eddy  diffusivily. 
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T = equilibrium  temperature 
k = Coriolis  parameter, 
f(AT)  = equation  of  state  and 

( ) = refers  to  dimensional  quantity. 

The  equilibrium  temperature  is  defined  as  the 
temperature  at  the  surface  for  which  there  is  no  heat 
transfer. 

The  conservative  form  of  the  convective  terms 
is  used  as  this  has  been  found  by  Arakawa^"^  to  be 
advantageous  for  numerical  computations.  The  density 
is  taken  as  only  a function  of  temperature.  The  energy 
equation  is  nondimensionalized  in  terms  of  temperature 
differences.  The  effect  of  round-off  error  will  be 
less  in  the  evaluation  of  the  derivatives  if  the  dif- 
ferences are  used. 

The  rigid-lid  condition  is  difficult  to  apply 
in  a numerical  solution  of  the  above  system  of  equa- 
tions. To  alleviate  this  difficulty,  an  additional 
equation,  a Poisson  equation  for  the  pressure,  which 
contains  the  rigid-lid  condition,  can  be  derived. 

This  is  accomplished  by  taking  the  divergence  of  the 
vertically  integrated  horizontal  momentum  equations 
and  using  the  vertically  integrated  continuity  and 
hydrostatic  pressure  equations.  The  Poisson  equation 
is : 
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The  term  P is  the  integration  constant  resulting 
from  the  vertical  integration  of  the  hydrostatic 
pressure  equation  and  is  the  surface  pressure,  i.e., 
the  pressure  at  the  surface  z = 0 . It  is  a function 
of  hoth  x and  y . This  surface  pressure  can  be 
interpreted  in  terms  of  the  equivalent  height  of 
water  above  or  below  the  surface  z = 0 required 
to  provide  the  prescribed  pressure.  In  this  way, 
surface  displacements  of  water,  neglecting  gravity 
waves,  can  be  compared  between  this  rigid-lid  model 
and  the  equivalent  free  surface  model. 

The  vertical  velocity  at  the  surface  z = 0 

has  not  been  set  to  zero  in  the  right-hand-side  of 

the  Poisson  equation  for  the  surface  pressure.  This 

r 32 1 

is  because  a corrective  procedure  is  used  in  the 
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numerical  solution  to  eliminate  accumulative  error 
in  the  satisfaction  of  the  continuity  equation. 

The  following  boundary  conditions  are  used  with 
the  above  system  of  equations: 


River  outflow 
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Bottom 


Surface 
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9v 

3n 


0 or 


v = 


f 


2 


9AT 

9n 


or  AT  = f. 


Pressure  Conditions 


9p 

= integrated  x or  y momentum  equation, 
specify  pressure  level  at  one  point. 

The  functional  forms  g^  , g^  , and  are 

the  specified  velocity  and  temperature  profiles  across 

the  river  outfall.  The  bottom  and  shore  are  taken  as 

no-slip,  impermeable,  insulated  surfaces.  A heat 

transfer  condition  proportional  to  the  temperature 

[ 32] 

difference  and  a wind-dependent  stress  are  imposed 
at  the  surface.  The  normal  derivative  pressure  boundary 
conditions  are  derived  from  the  appropriate  vertically 
integrated  momentum  equation.  The  pressure  must  be 
specified  at  one  point  to  make  its  solution  unique. 

The  boundary  conditions  applied  at  the  outer  x and 
y boundaries  are  either  that  the  normal  derivatives 
of  the  velocity  and  temperature  are  zero,  or  that  the 
velocity  and  temperature  are  specified. 
lU6.  The  general  arrangement  of  variables  in  the  grid  system  for 
the  plume  model  is  almost  identical  with  that  used  previously  by  Paul 
and  Lick.  ~ The  horizontal  velocities  are  defined  at  integral  nodal 
points,  the  vertical  velocity  is  defined  at  half-integral  nodal  points, 
the  temperature  is  defined  at  half-integral  nodal  points  in  the  hori- 
zontal and  integral  nodal  points  in  the  vertical,  and  the  surface  pres- 
sure is  defined  at  half-integral  nodal  points  in  the  horizontal.  Fig- 
ure 29  indicates  typical  horizontal  and  vertical  grid  sections  and 
shows  the  relative  positions  of  the  various  variables.  For  the  deri- 
vation of  the  finite  difference  equations,  variables  are  sometimes 
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Figure  29.  Arrangement  of  variables  in  grid  sections 


required  at  points  where  they  are  not  defined.  In  these  circumstances, 
the  undefined  quantity  is  taken  as  the  simple  average  of  the  neighboring 
values . 

ll*7.  The  Euler  explicit  time  scheme  is  used  exclusively  in  the 
present  model  in  which  the  time  derivative  is  written  as  a simple  for- 
ward time  difference  and  the  rest  of  the  equation  is  evaluated  with  the 
previously  calculated  values. 

lU8.  The  following  scheme  is  used  for  the  solution  of  the  differ- 
ence equations  for  the  plume  model: 

a.  It  is  assumed  that  values  from  the  previous  time  step  are 
available . 

b.  The  surface  pressure  is  calculated  with  the  right  side 
of  the  equation  evaluated  from  the  previous  time  step 
values . 

£.  The  temperature  is  calculated  by  the  explicit  time 
scheme . 

d_.  The  density  is  calculated  from  the  equation  of  state. 

e_.  The  horizontal  velocities  are  calculated  by  the  explicit 
time  scheme. 

f\  The  vertical  velocity  is  calculated  by  vertically  inte- 
grating the  continuity  equation  from  the  bottom. 

£.  The  present  time  step  is  now  complete. 

At  each  time  step  the  Poisson  equation  in  the  surface  pressure  has  been 
solved.  This  is  solved  by  the  alternating-direction-implicit  (ADI) 
method.33,31* 

11*9.  After  the  temperatures  are  calculated  by  the  explicit  time 
scheme,  they  are  checked  to  see  if  static  stability  is  satisfied,  i.e., 
if  the  temperatures  decrease  monotonically  downward  (assuming  that 
density  increases  with  decreasing  temperature).  When  a static  insta- 
bility is  encountered,  an  infinite  mixing  procedure  is  used.  This 
procedure  merely  averages  the  temperature  over  an  unstable  region. 

190.  The  thermal  plume  model  has  been  applied  for  the  Point  Beach 
Power  Plant  thermal  outfall  and  for  the  Cuyahoga  River  entering  Lake 
Erie.  In  the  vicinity  of  the  inlet,  the  grid  spacing  is  constant,  then 
increases  with  distance  from  the  inlet,  and  finally  becomes  constant  in 
each  direction.  Such  a grid  system  is  used  to  pick  up  details  of  the 
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flow  near  the  inlet  where  quantities  change  rapidly  and  to  avoid  an 
excessive  number  of  points  away  from  the  inlet  where  quantities  are 
changing  much  more  slowly.  Constant  spacing  is  used  in  the  vertical 
direction . 

151.  The  vertical  eddy  coefficient  is  taken  as  dependent  on  the 
local  vertical  temperature  gradient  (3T/3z).  The  actual  form  used  is 
similar  to  that  used  in  an  application  of  the  numerical  model  to  the 
wind-driven  and  thermally  driven  circulation  in  a small  pond.^'’ 


where  a and  6 are  constants  dependent  on  the  local  conditions  of 
the  problem.  The  constant  a is  chosen  so  that  in  the  absence  of 
temperature  gradients,  the  eddy  coefficient  is  equal  to  that  which  would 
be  used  in  a constant  eddy  coefficient  model. 

152.  The  stress  acting  on  the  water  surface  due  to  the  wind  is 

2I4 

calculated  by  the  formulae  developed  by  Wilson.  These  formulae  have 

been  used  successfully  in  numerical  calculations  of  wind-driven  circu- 

22  23  35 

lations  in  lakes  and  small  ponds.  ’ ’ 

153.  The  results  obtained  from  the  plume  model  for  these  two 
indicated  applications  appear  reasonable;  however,  no  comparison  with 
prototype  data  is  available.  The  complete  lake  and  plume  model  has  not 
as  yet  been  applied  to  any  actual  physical  geometries. 


I 
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PART  V:  RECOMMENDED  HYDRODYNAMIC  MODELS 


Storm  Surge 

15l*.  From  a review  of  available  literature,  accurate  storm  surge 
surface  elevations  in  lakes  apparently  can  be  obtained  from  the  two- 
dimensional  depth-averaged  models  if  a free  surface  boundary  condition 
is  used.  If  hydrodynamic  quantities  other  than  surface  elevations  are 
desired,  however,  the  two-dimensional  storm  surge  model  is  of  little 
use  and  three-dimensional  models  must  be  used.  The  changes  in  circu- 
lation patterns  produced  by  structures  can  be  obtained  only  from  three- 
dimensional  models. 

155.  The  three-dimensional  models  are  usually  more  expensive  to 
operate  than  the  two-dimensional  models.  This  expense  suggests  the 
use  of  two-dimensional  models  in  the  preliminary  storm  surge  calcu- 
lations followed  by  application  of  the  three-dimensional  models  for 
final  detailed  analyses. 

156.  Several  satisfactory  two-dimensional  storm  surge  models  with 

a free  surface  condition  are  available.  The  three-dimensional  storm 

2 3 

surge  model  developed  by  Haq,  Lick,  and  Sheng  appears  the  most  promis 
ing  of  a limited  number  of  three-dimensional  models  and  should  be  avail 
able  in  a matter  of  months  with  initial  application  to  Lake  Erie  for 
storm  Agnes . 


Wind-Driven  Circulation 


157.  The  two-dimensional  models  provide  very  lit.tli  u:'"fu!  inf  • - 

mation  relative  to  lake  circulation.  The  velocity  pre-1 1 -t  i i,  1 • n-  i 
from  the  two-dimensional  models  are  interesting  only  in  l r . , j - 
tative  manner.  The  detailed  hydrodynamic  variable:-  r»-  ; . r • i ' 1-  • • 

mine  the  complex  velocity  structure  in  a lake  can  he  t • -i  i n>  1 r . • ■ 

a three-dimensional  model. 

158.  The  three-dimensional  numerical  model:-  {recently  • » v • 1 : 

and  in  the  development  stage  by  Dr.  Wilbert  Lick  an  i 1 r 
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Case  Western  Reserve  University  appear  to  be  the  most  promising  models 
for  wind-driven  circulation  in  freshwater  lakes.  These  models  are 
particularly  appropriate  for  Lake  Erie  applications  since  the  models 
have  previously  been  applied  to  Lake  Erie  and  much  lake  information  is 
already  available  in  a convenient  format  for  input  into  the  model. 

159*  A three-dimensional  steady-state  circulation  model  is  pres- 

22 

ently  operational  at  Case  Western  Reserve  University.  This  model  has 
been  applied  to  Lake  Erie  with  some  degree  of  verification  based  on 
prototype  data.  A three-dimensional  stratified  circulation  model  is 
presently  under  development  and  should  be  operational  in  the  near 
future. 


PART  VI:  MODEL  VERIFICATION 


General  Verification  Procedures 


160.  A hydrodynamic  mathematical  model  is  a discretization  of  the 
actual  physical  system.  As  such,  the  mathematical  model  is  somewhat 
limited  in  the  resolution  of  actual  physical  details.  In  a mathematical 
model  based  upon  finite  difference  techniques,  any  physical  detail 
smaller  than  the  grid  size  cannot  be  directly  represented.  In  a similar 
manner,  the  results  of  the  mathematical  model  cannot  resolve  small-scale 
local  circulation  patterns  or  eddies  unless  the  condition  extends  over 
an  area  covering  several  grid  points.  In  addition,  since  the  depths  are 
indicated  in  a discrete  manner,  the  numerical  model  tends  to  smooth  or 
average  the  bottom  topography  and  eliminate  local  undulations.  At  the 
same  time,  changes  in  the  boundary  configuration  are  made  in  a step 
manner  rather  than  a gradual  change. 

161.  The  mathematical  modeler  must  therefore  select  a grid  size 
and  orientation  of  the  coordinate  system  that  will  allow  the  desired 
phenomena  to  be  resolved.  He  must  also  use  judgment  and  experience  in 
the  selection  of  system  boundaries  and  the  representation  of  the  topog- 
raphy. These  details  represent  part  of  the  "art"  as  opposed  to 
"science"  of  mathematical  modeling  and  must  be  accomplished  prior  to 
actually  "running"  the  model.  Unfortunately,  although  very  important, 
these  details  are  some  of  the  lesser  problems  associated  with  mathe- 
matical models. 

162.  The  mathematical  model  formulation  contains  parameters  that 
are  not  well  defined  quantities.  Bottom  friction,  the  horizontal  eddy 
viscosity,  the  vertical  eddy  viscosity  coefficient,  and  the  wind  speed- 
surface  stress  relationship  are  examples  of  model  input  that  are  very 
much  dependent  upon  the  specific  model  application.  The  bottom  friction 
is  basically  a function  of  the  bottom  surface,  i.e.  sand,  mud,  rubble, 
grass,  etc.  The  eddy  viscosity  coefficients  are  functions  of  the  degree 
of  turbulence,  the  density  variations,  the  wind  condition,  the  local 
ship  traffic,  etc.  The  wind  speed-surface  stress  relationship  presents 
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problem?  primarily  due  to  difficulties  associated  with  extrapolating 
known  wind  conditions  at  some  fixed  land  site  to  obtain  a wind  speed 
at  some  location  in  the  lake. 

163.  The  nature  of  the  parameters  in  the  mathematical  model  are 
such  that  model  verification  is  required  prior  to  using  the  model  as  a 
predictive  tool.  Order  of  magnitude  values  of  the  parameters  are  known; 
however,  these  values  are  not  sufficient  for  the  model  applications. 

The  model  parameters  must  be  adjusted  in  such  a manner  that  the  numeri- 
cal model  results  agree  with  existing  prototype  conditions.  Once  agree- 
ment with  existing  prototype  condition  has  been  obtained,  the  values  of 
the  parameters  required  to  obtain  the  agreement  are  considered  to  be 
the  appropriate  values  of  the  parameters.  These  values  of  the  parame- 
ters are  then  used  when  the  model  is  applied  -to  .predict  changes  result- 
ing from  some  proposed  modification  of  the  existing  condition. 

l6^.  Sufficient  prototype  data  for  existing  conditions  must  be 
known  to  allow  the  mathematical  model  to  be  verified.  The  degree  of 
verification,  and  hence  confidence,  in  the  model  results  depend  upon 
the  amount  of  prototype  data  available  for  verification  of  the  model. 

For  complete  model  verification,  prototype  data  should  exist  for  all 
hydrodynamic  variables  under  investigation,  and  the  data  should  be 
available  at  a number  of  locations  extending  over  the  area  to  which  the 
model  will  fee  applied.  If  less  than  this  amount  of  prototype  data  is 
available  for  model  verification,  a corresponding  decrease  in  confidence 
in  the  model  results  must  be  accepted. 

Verification  Status  of  Recommended  Models 


Storm  surge 

165.  The  two-dimensional  storm  surge  models  are  the  most  exten- 
sively verified  of  the  recommended  models.  While  there  are  some  differ- 
ences between  the  various  models,  the  basic  formulation  procedure  is 
the  same  and  similar  results  are  obtained  from  all  the  models.  The 
surface  elevations  predicted  by  the  two-dimensional  models  have  been 
compared  with  measured  values  for  several  actual  storm  systems.  For 
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example,  Platzman  has  investigated  several  extreme  storm  cur Re  con- 

20 

ditions  on  Lake  Erie;  Paskausky  has  investigated  the  Lake  Erie  storm 

l8 

surge  associated  with  storm  Agnes;  and  Reid  and  Bodine  have  investi- 
gated storm  surge  in  Galveston  Bay  during  hurricanes  Carla  (1961)  and 
Cindy  (1963).  The  surface  elevations  predicted  by  the  two-dimensional 
models  display  a reasonable  agreement  with  measured  tidal  elevations. 

166.  The  recommended  three-dimensional  storm  surge  model  has  been 

23 

partially  verified  by  application  to  Lake  Erie  during  storm  Agnes. 

The  calculated  and  measured  values  for  the  surface  displacement  at 
Cleveland  and  at  Port  Stanley  were  found  to  be  in  good  agreement . There 
has  been  no  known  effort  to  verify  the  calculated  velocities,  due  pri- 
marily to  a lack  of  reliable  prototype  velocity  data. 

Wind-driven  circulation 

167.  The  recommended  three-dimensional,  steady-state,  constant- 
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density  model  of  Gedney  and  Lick  has  been  verified  to  a limited 
extent.  The  model  has  been  applied  to  Lake  Erie  for  a constant  wind 
condition  and  the  calculated  velocities  compared  with  measured  veloci- 
ties observed  whenever  the  wind  was  fairly  steady  for  two  or  more  days. 
The  results  show  that  the  velocities  vary  greatly  from  position  to 
position  and  depend  strongly  on  the  bottom  topography  and  boundary 
geometry.  The  calculated  velocities  are  compared  with  a limited  number 
of  continuous  current  meter  measurements  at  10  and  15  m below  the  sur- 
face. Good  quantitative  agreement  between  the  meter  measurements  and 
calculations  were  found  to  exist  whenever  the  wind  was  fairly  steady  for 
two  or  more  days.  Qualitative  agreement  with  lakebed  drifter  and  sur- 
face drift  measurements  has  also  been  observed, 

168.  The  three-dimensional,  variable-density  model,  which  is 
recommended,  is  not  yet  completely  operational;  and  no  model  verifi- 
cation efforts  have  been  conducted.  This  model  is,  however,  an  exten- 
sion of  a thermal  plume  model  to  which  some  verification  efforts  have 
been  directed. 

169*  The  recommended  models  have  not  been  verified  to  the  extent 
needed  for  detailed  quantitative  results.  This  situation  is  particu- 
larly true  for  the  velocity  fields  predicted  by  the  three-dimensional 
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models.  Additional  Lake  Erie  prototype  data,  especially  for  the  region 
near  Cleveland,  are  required  to  verify  the  models  and  to  enhance  their 
operational  mode  as  a predictive  tool. 
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PART  VII : RECOMMENDED  PROTOTYPE  DATA  ACQUISITION  SYSTEM 


Prototype  Data  Requirements 

170.  A review  of  available  data  by  WES  as  a part  of  the  Lake  Erie 
International  Jetport  Model  Feasibility  Investigation  indicated  that 
sufficient  prototype  data  near  Cleveland  were  not  available  to  verify 
either  a physical  or  analytical  lake  mass  circulation  model.  Other 
areas  where  insufficient  data  exist  include  wave  regime  definition, 
littoral  transport,  and  harbor  flushing.  Limited  current  data  are 
available  at  four  stations  near  Cleveland  for  several  months  during  the 
summer  of  1965,  from  a 196U-65  Environmental  Protection  Agency  (formerly 
the  Federal  Water  Pollution  Control  Agency)  current  metering  program. 

The  data  from  these  stations  are  not  continuous  over  the  recording 
period  and,  for  three  of  the  stations,  are  available  for  one  depth  only. 

171.  A prototype  data  acquisition  and  analysis  program  is  proposed 

by  WES  to  obtain  sufficient  prototype  data  for  verification  of  analyti- 
cal and  physical  circulation  models  of  the  effects  of  a jetport  on  the 
hydrodynamics  of  Lake  Erie  near  Cleveland,  Ohio.  This  program  has  two 
modes:  (a)  time-series  data  in  the  lake  proper  and  (b)  synoptic  data 

near  and  in  the  commerical  harbor  and  river  navigation  channel.  For  the 
time-series  mode,  water  velocity,  water  temperature,  mean  water  level, 
wave  height,  wind  velocity,  and  air  temperature  will  he  monitored;  and 
in  the  synoptic  survey,  water  velocity,  water  temperature,  and  dye 
concentration  will  be  monitored. 

Proposed  Time-Series  Data  Acquisition  System 

172.  In  the  time-series  mode  of  data  acquisition,  the  hydrological 
parameters  to  be  recorded  are  the  x-y  components  of  water  and  wind  ve- 
locities, water  and  air  temperature,  wave  height,  water  elevation,  and 
barometric  pressure.  A data  acquisition  system  to  collect  these  data 
can  be  subdivided  into  the  transducer  conditioner,  transmitter,  re- 
ceiver, conditioner  for  recording,  and  recorder  subsystems. 
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173.  Four  instrument  towers  are  proposed  to  be  installed  in  Lake 
Erie  near  Cleveland.  The  transducer,  conditioner,  and  transmitter  sub- 
systems are  to  be  mounted  on  these  towers  and  shall  be  capable  of  oper- 
ating at  rated  specifications  in  the  ambient  environmental  conditions. 
Location  of  the  towers  and  the  instrumentation  to  be  attached  to  each 
tower  are  shown  in  Figures  30  and  31.  The  towers  will  be  similar  in 
configuration  and  bracing  to  triangular  communication  towers  mounted  on 
land.  They  will  be  80  ft  in  height  and  installed  in  a water  depth  of 
about  60  ft.  Two  sets  of  guy  wires  will  be  included  as  shown  in  Fig- 
ure 32,  The  power  requirements  shall  be  kept  at  a minimum  to  ensure 
the  longest  battery  life  possible.  The  analog  output  of  the  velocities 
and  temperature  transducers  shall  be  digitized  and  conditioned  for  radio 
transmission . 

17^.  The  receiver,  conditioner  for  recording,  and  recording  sub- 
systems shall  be  located  in  an  office  area  where  air  conditioning, 
regulated  power,  and  general  laboratory  conditions  are  available.  These 
subsystems  need  meet  only  laboratory  type  specifications. 

Specifications  for  transducers 

175.  The  transducer  subsystem  will  require  various  sensors  that 
are  specialized  to  collect  the  required  data  at  specified  sampling  rates. 
Such  instruments  are  highly  specialized  and  will  be  discussed  below  as 
to  the  required  data  types. 

176.  Water  velocity.  A magnetic  current  meter  capable  of  measur- 
ing the  x-y  components  of  the  water  velocity  shall  be  used.  This  type 
meter  was  selected  because  of  its  small  size,  accuracy,  and  no  moving 
parts,  which  would  require  constant  attention.  The  meter  range  is  0 to 
+_  10  fps  +_  0.07  fps  or  2 percent  of  reading.  Increased  accuracy  to 

+ 0.02  fps  can  be  obtained  at  additional  cost.  The  resolution  shall -be 
0.03//T  fps  where  T is  the  time  constant.  The  standard  time  constant 
is  1 sec  with  optional  values  from  0.2  sec  to  20  sec.  There  shall  be 
a recorder  output  of  + 5 v full  scale.  The  current  meter  sensor,  which 
will  be  mounted  vertically,  is  approximately  10-1/2  in.  long  and  1 in. 
in  diameter. 

177.  Wind  velocity.  The  wind  velocity  shall  be  measured  with  two 
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Location  of  time-series  survey  stations 


a.  TOWER  LOCATIONS 


TOWERS  A,  C,  S D TOWER  B 


b.  INSTRUMENTATION  LOCATIONS 

Figure  31.  Location  of  towers  and  instruments  for  time-series  survey 


i 


Figure  32.  Schematic  of  instrumentation  towers 
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perpendicularly  mounted  propeller-type  anemometer:-,  that  respond  only  to 
that  component  of  wind  parallel  with  its  axis  of  rotation.  When  the 
wind  is  exactly  perpendicular  to  the  axis  of  the  propeller,  it  will  stop 
rotating.  The  propeller  respionds  as  a function  of  its  orientation  to 
the  wind;  this  response  very  closely  approximates  the  cosine  law.  Each 
sensor  measures  both  forward  and  reverse  air  movement,  providing  a 
signal  whose  magnitude  is  proportional  to  propeller  speed  and  whose 
polarity  indicates  direction.  The  range  shall  be  90  mph  head  on  and 
TO  mph  at  all  angles.  The  threshold  shall  be  0.3-0. 5 mph.  The  pro- 
peller shall  have  0.96  revolutions  per  foot  of  air  movement.  The 
d = c tachometer  output  shall  be  + 500  mv  DC  at  +1800  rpm.  u 

178.  Water  elevation.  The  water  elevation  transducer  shall  con- 
sist of  a mechanically  driven  encoder,  a float,  tape,  pulley,  and 
counterweight.  The  float  shall  operate  in  a cylinder  that  is  vented 
to  the  lake  through  a small  hole  near  the  bottom  of  the  cylinder.  The 
water  level  in  the  cylinder  is  the  mean  water  elevation  of  the  lake  and 
will  not  vary  with  wave  action.  As  the  water  level  changes,  the  float 
system  rotates  the  encoder  input  shaft.  The  encoder  shall  have  binary 
coded  decimal  output  with  a resolution  and  accuracy  of  one  part  in 
10,000.  Full  range  shall  be  99-99  ft. 

179-  Water  temperature.  The  temperature  sensor  shall  be  capable 
of  measuring  temperature  over  a range  of  0°C  to  bo°C  to  an  accuracy  Oi 
+0.15°C.  The  output  of  the  device  shall  be  no  less  than  1 v full 
scale. 

130.  Air  temperature.  The  air  temperature  sensor  shall  be  capable 
of  measuring  temperature  over  a range  of  -10°C  to  30°C  to  an  accuracy 
of  +0.25°C.  The  sensor  shall  be  mounted  in  a structure  to  shield  it 
from  direct  sun.  The  output  shall  be  5 v full  scale. 

l8l.  Wave  height.  Wave  heights  will  be  measured  with  two  trans- 
ducers, one  at  a shallow  depth  for  small  waves  and  one  at  a greater 
depth  for  large  waves.  Each  wuve  height  transducer  will  be  an  absolute 
pressure  transducer.  This  transducer  will  have  no  moving  parts  exposed 
to  the  lake  water.  The  pressure  transducer  will  be  encased  in  a con- 
tainer 6.5  in.  0D  and  15  in.  long.  The  pressure  port  will  be  in  the 
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top  of  this  container,  and  the  case  will  be  made  of  plastic  or  some  non- 
corrosive,  nonmetallic  material.  The  shallow-depth  pressure  transducer 
range  will  be  58  ft  absolute,  from  which  the  barometric  pressure  must  be 
subtracted,  leaving  about  25  ft  for  wave  variations.  The  greater-depth 
transducer  range  will  be  116  ft  absolute,  with  about  80  ft  for  wave 
variations.  The  operating  temperature  will  be  35°F  to  T5°F.  The  output 
will  be  0.28  to  7-5  MHz  with  linearity  and  hysteresis  being  0.05  percent 
of  full  scale  from  best  straight  line  between  15  psia  and  full  scale. 

The  repeatability  will  be  0.02  percent  full  scale  with  a long-term  accu- 
racy of  0.25  percent.  The  frequency  analog  output  will  be  connected  to 
an  above- the-sur face  unit  that  will  convert  the  frequency  to  a 12-bit 
natural  binary  output. 

182.  Sampling  rates.  The  sampling  rates  of  each  instrument  will 
vary  depending  upon  the  phenomena  to  be  sampled.  The  anticipated  sam- 
pling rates  are  as  follow: 
a_.  Temperature. 


(1) 

Air: 

instantaneous  value  every  20  min. 

(2) 

Water : 

: instantaneous  value  every  20  min. 

b. 

Velocity. 

(1) 

Air : 

20-min  average  of  1-sec  digital  valves. 

(2) 

Water : 

: 20-min  average  of  1-sec  digital  valves 

c. 

Wave 

s • 

(1) 

Deep : 

four  per  second  for  20  min  every  3 hr. 

(2)  Shallow:  four  per  second  for  20  min  every  3 hr. 

ci.  Lake  level.  Instantaneous  value  every  20  min. 

183.  Analog-to-digital  converter  (ADC).  The  ADC  converts  the 

analog  voltage  outputs  of  the  various  transducers  to  a 2's  complement 

binary  code.  A l6-channel  analog  multiplexer  will  be  an  integral  part 

of  the  ADC  unit.  The  analog  input  levels  are  +5  v.  The  input  im- 

9 

pedance  will  be  10  ohms.  The  transfer  characteristics  will  have  an 
accuracy  and  linearity  of  0.023  percent  full  scale  + least  significant 
bit.  The  temperature  coefficient  will  be  no  greater  than  O.OUI4  per- 
cent/°C.  The  acquisition  time  will  be  no  greater  than  150  psec.  The 
system  throughput  rate  will  be  no  greater  than  I450  psec.  The  ADC  unit 


will  operate  from  a 12-v  storage  battery.  The  binary  output  will  be 
positive  true  logic  and  be  transistor-true-logic  compatible. 

18U . Conditioning  and  transmission.  The  output  of  the  ADC  will  be 
connected  to  shift  registers  where  it  will  be  stored  until  ready  for 
radio  telemetry.  The  average  values  of  velocity  could  be  made  at  the 
transducer  location  and  then  conditioned  and  transmitted  by  radio  te- 
lemetry to  the  central  station.  Alternatively,  the  average  velocity 
values  could  be  made  after  transmission,  if  determined  to  be  more  feasi- 
ble and  economical.  The  radio  telemetry  network  shall  have  the  capa- 
bility of  multiplexing  all  data  channels,  either  parallel  or  serial. 
Block  diagrams  of  these  systems  are  shown  in  Figures  33  and  3^. 

185.  Central  station.  The  telemetry  signals  are  received  and  con- 
verted to  a digital  format  compatible  with  a 1/2-in.  magnetic  tape  re- 
corder. These  magnetic  tapes  will,  in  turn,  be  compatible  with  a large 
General  Electric  computer  at  WES  and  will  be  shipped  to  Vicksburg, 
Mississippi,  at  regular  intervals  for  analysis  and  data  reduction.  All 
of  the  parameters  shall  be  recorded  on  the  magnetic  tape.  There  shall 
be  a digital  clock  which  shall  keep  time  in  day  of  the  year,  hour, 
minute,  and  second.  This  clock  shall  have  a binary  coded  decimal  out- 
put compatible  with  the  magnetic  tape  recorder.  The  clock  shall  keep 
time  for  the  entire  system.  From  its  output,  timed  sync  signals  shall 
be  transmitted  to  the  remote  stations.  These  sync  signals  shall  time 
the  acquisition  and  transmission  of  data. 

186.  The  velocity,  temperature,  water  elevation,  barometric 
pressure,  and  time  data  shall  be  printed  out  on  paper  for  a "quick  look" 
to  determine  if  these  channels  are  working.  The  water  wave  data  will 

be  recorded  on  a strip  chart  recorder.  The  strip  chart  recorder  shall 
be  turned  on  and  off  with  the  systems  clock.  A block  diagram  of  the 
central  station  is  shown  in  Figure  35. 

Summary 

187.  The  transducers  have  been  given  more  time  and  consideration 
because  it  is  felt  that  methods  of  detection  and  sensors  are  of  prime 
importance.  The  technique  and  hardware  for  the  conversion,  handling, 
transmission,  reception,  and  recording  of  these  data  can  be  designed 
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LOCATION  B 
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Figure  35.  Block  diagram  base  station 


fully  after  a prototype  study  has  been  authorized.  The  block  diagrams 
(Figures  33-35)  depicting  the  initial  concept  of  the  acquisition  system 
and  the  receiving  and  recording  system  at  the  base  station  are  illustra- 
tive in  purpose,  are  preliminary  in  nature,  and  do  not  constitute  the 
final  design  of  the  systems. 

Proposed  Synoptic  Survey 

188.  The  parameters  to  be  monitored  in  each  of  the  two  synoptic 
surveys  are  current  velocity,  water  temperature,  and  dye  concentration. 
Observations  of  these  parameters  will  be  taken  at  2-hr  intervals  over 
a 12-hr  period  for  8 days.  These  data  will  be  monitored  at  three  depths 
for  32  stations.  Location  of  these  stations  is  indicated  in  Figure  36. 
The  stations  will  be  marked  by  buoys  anchored  to  the  lake  bottom.  The 
surveys  will  be  made  during  a typical  spring  or  fall  condition  and 
during  a typical  summer  stratified  condition. 


Figure  36.  Location  of  synoptic  survey  stations 
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189.  Pye  concentrations  during  the  survey  will  be  obtained  after 
point  release  of  a known  quantity  of  dye  in  the  Cuyahoga  River  near  its 
outflow  into  the  harbor.  The  concentrations  in  water  samples  from  each 
station  will  be  determined  using  a fluorometer.  Current  velocity  will 
be  measured  using  low-threshold  electromagnetic  current  meters  (required 
for  the  low  magnitudes  of  currents  typical  of  the  lake  circulation). 

190.  The  prototype  data  monitored  during  the  survey  will  be  tran- 
scribed from  field  observations  to  a format  compatible  with  the  WES 
computer  prior  to  analysis  of  the  data.  Analyses  of  these  data  basi- 
cally consist  of  determining  the  flushing  rate  and/or  time  required  for 
conservative  constitutents  to  move  out  of  the  present  harbor,  estimating 
local  horizontal  and  vertical  eddy  diffusivity  coefficients,  and  obtain- 
ing statistical  information  on  the  movement  of  the  outflow  plume  of  the 
Cuyahoga  River. 

191 • As  a part  of  the  survey,  surface  temperature  patterns  and,  if 
present,  the  location  of  the  Cuyahoga  River  sediment  plume,  could  be 
obtained  by  overflights  conducted  by  the  NASA  Lewis  Research  Center. 
Initial  inquiries  by  WES  have  been  conducted  with  NASA  to  define  WES's 
requirements  and  NASA's  support  capabilities  for  this  survey.  Initial 
efforts  indicated  an  existing  potential  for  cooperative  efforts,  but 
further  study  definition  and  coordination  will  be  required  in  advance  of 
any  field  study. 


PART  V I T 1 : CONCLUSIONS  AND  RECOMMENDATIONS 


192.  Current  state-of-the-art  numerical  models  and  numerical 
models  that  may  become  available  in  the  time  frame  necessary  for  appli- 
cation to  the  Lake  Erie  Jetport  study  can  provide  much  useful  hydrody- 
namic information.  The  presently  available  seiche  models  and  two- 
dimensional  wind-driven  circulation  models  have  been  operational  for 
several  years  and  within  their  recognized  limitations  appear  to  yield 
valid  information.  A three-dimensional,  constant-density,  steady-flow 

model  for  Lake  Erie  was  reported  in  the  literature  in  1972  by  Gedney 

22 

and  Lick.  Dr.  Wilbert  Lick  and  associates  at  Case  Western  Reserve 

23 

University  are  working  on  three-dimensional  storm  surge  models  and 
three-dimensional  variable-density  models. The  three-dimensional 
storm  surge  models  should  be  available  in  an  operational  form  in  a 
matter  of  months.  The  time  schedule  on  the  three-dimensional,  variable- 
density  model  is  less  well  defined  due  to  the  complexity  of  the  model; 
however,  it  can  hopefully  be  made  available  within  the  time  frame  when 
application  to  the  jetport  study  will  be  possible. 

193.  In  general,  the  numerical  models  have  not  been  verified  to 
the  extent  desirable.  If  they  are  applied  to  the  jetport  study,  addi- 
tional model  verification  will  be  necessary  for  which  sufficient  proto- 
type data  are  not  available. 

19^.  The  following  recommendations  are  the  result  of  the  numerical 
model  feasibility  study: 

a.  Obtain  for  the  jetport  study  operational  mathematical 
models  which  have  been  recommended  earlier  in  this 
report . 

b.  Insure  the  expedient  development  of  needed  three- 
dimensional  models,  which  are  presently  in  a development 
state,  by  contracting  with  Case  Western  Reserve  Uni- 
versity to  provide  these  models  with  preliminary  appli- 
cation to  the  jetport  study. 

c_.  Provide  the  necessary  prototype  data  for  proper  verifi- 
cation of  numerical  models  (and/or  future  physical 
models)  by  initiating  the  prototype  data  acquisition 
program  outlined  earlier  in  this  report. 

d.  Estimate  the  effects  of  a proposed  offshore  Jetport  on 


storin  surge  and  lake  circulation  in  Lake  Erie,  particu- 
larly in  the  vicinity  of  Cleveland,  Ohio,  by  conducting 
the  following  recommended  numerical  studies: 

(1)  Investigate  the  near-field  and  far-field  effects  of 
a jetport  island  on  the  steady-state  wind-driven 
circulation  in  Lake  Erie  for  well  mixed  (nonstrati- 
fied)  lake  conditions.  Circulation  patterns  with 
and  without  a jetport  should  be  defined  for  average 
and  extreme  winter  wind  speeds  and  several  wind 
directions. 

(2)  Investigate  the  effects  of  a jetport  island  on  the 
storm  surge  along  the  shoreline  from  Fairport  to 
Lorain,  Ohio.  Time-dependent  wind  fields  for  at 
least  two  severe  storms  on  Lake  Erie  should  be  used 
to  estimate  surge  elevation  and  velocity  fields  for 
well  mixed  lake  conditions. 

(3)  Investigate  the  effects  of  a jetport  island  on  the 
stratified  lake  circulation  (epilimnion  and  hypo- 
limnion)  around  the  jetport  and  in  the  Cleveland 
Harbor  with  typical  Cuyahoga  River  outflows  and 
steady-state,  summer  wind  fields. 
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Table  1 


Relative  Amplitude  Comparison 


Location 

Observed  Relative 
Amplitude* 

Computed  Relative 
Amplitude 

Buffalo 

-0.70 

-0.72 

Port  Colborne 

-0. 66 

-0.69 

Dunkirk 

-0.66 

-0 . 66 

Port  Dover 

-o.6o 

-0.63 

Erie 

-0.50 

-0.1*8 

Port  Stanley 

-0.19 

-0.09 

Fairport 

0.13 

0.09 

Erieau 

0.27 

0.21 

Cl eveland 

0.25 

0.27 

Huron 

0.1*9 

0.54 

Port  Clinton 

0.52 

0.86 

Toledo 

1.00 

1.00 

Monroe 

0.95 

0.99 

\ 

* Based  on  results  of  spectral  analysis 
variations  (Reference  9)- 

of  observed  lake  level 

Periods  for 

Table  2 

First  Five  Modes 

of  Oscillation 

for 

Different  Jetport  Configurations 

Period,  hr 

Mode  of  Oscillation 

Configuration* 

1 

2 

3 

4 

5 

No  Jetport — Observed 

lU . 38 

9.14 

5.93 

4.15 

— 

No  Jetport — Calculated 

ll*  .1*3 

9.22 

6.01 

4.31 

3.87 

A 

ll* . 50 

9.28 

6.06 

4.32 

3.87 

B 

11+.51 

9.28 

6.05 

4.33 

3.87 

C 

11+.52 

9.28 

6.07 

4.33 

3.87 

D 

14.52 

9.27 

6.06 

4.33 

3.87 

* Letters  refer  to  configurations  shown  in  Figure  5* 
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APPENDIX  A:  NOTATION 


A 

[A] 


A 

A, 


o 

B 

C 

[C] 

c* 

cd 

c 

p 

D 


D ,D  ,D 
x y z 


V 

f 

f’  (T) 
f 

o 

g 

G 

w 

h 


h(x,y ) 


Surface  area  of  Lake  Erie 

Square  matrix  whose  order  is  equal  to  the  number  N of  the 
interior  grid  points 

A divided  into  N subregions  (elements) 

Horizontal  eddy  diffusion  coefficient 
Vertical  eddy  diffusion  coefficient 

Horizontal  reference  width  of  river  discharging  the  plume 
Bottom  stress  coefficient 
Chezy  coefficient 

Vector  column  matrix  containing  the  nonhomogeneous  portion 
of  the  difference  equations 

Resistance  coefficient 

Drag  coefficient 

Specific  heat  at  constant  pressure 

Frictional  depth;  water  depth  in  Paskausky's  model 

Eddy  diffusion  coefficients  of  salinity  in  x , y , and 
z direction 

Horizontal  Ekman  number 

Vertical  Ekman  number 

Coriolis  parameter 

Equation  of  state 

Coriolis  parameter  evaluated  to  mean  latitude  of  Lake  Erie 

Acceleration  due  to  gravity 

Element  slope  matrix 

H/H  , dimensionless  depth 
n 

Water  depth 


A1 


he 


h 

0 
H 

[H] 

H 

m 

1 

-t 

1 


i,J  *k 

3 

k 


K 

[K] 


K 

v 

L 


m 

M 

[M] 

“t 

M 

x 

M 

y 

n 

n 

[N] 


Mean  depth  for  each  element 

Vertical  reference  depth  of  river  discharging  the  flume 

Characteristic  length  in  the  vertical  direction 

M x M square  matrix  of  coefficients 

Mean  depth  of  Lake  Erie 

Imaginary  number 

Unit  vector  in  x direction 

Indices  of  nodes  of  each  element 

Unit  vector  in  y direction 

Thermal  conductivity 

Coefficient  of  eddy  viscosity  or  diffusivity 
Element  of  stiffness  matrix 
Horizontal  diffusion  coefficient 
Vertical  diffusion  coefficient 

Characteristic  horizontal  scale;  characteristic  length  in 
the  horizontal  direction 

Nodal  point 

Total  number  of  node  points 
Lumped  mass  matrix 
Horizontal  mass  flux 
x-component  of  horizontal  mass  flux 
y-component  of  horizontal  mass  flux 
Unit  normal  to  boundary 
Manning  friction  factor 
Interpolation  or  shape  function 


A2 


1 


p 

p 

pk 

Q 

Q* 

Re 


R 

o 


S 

t 

R 

T 

T 


T 

s 


u 

u 

u' 


U 


V 


Value  of  interpolation  function  at  element  node  with 
index  i and  is  equal  to  [(x^y.  - x.y^)  + 


(yk  - 7^*  + 

(x.  - x )y]/2A  where  A is  the  area  of  the  element  k'  . 

,1  K 


Similar  expressions  obtained  by  the  cyclical  permutation  of 
i , j , and  k 


Pressure 


Volumetric  flow  in  the  x direction 

Layer-averaged  pressure  of  the  k'*'*1  layer 

Volumetric  flow  in  y direction 

Heat  source  and/or  sink 

Real  number 

Rossby  number 

Salinity 

Time 

Reference  time  value 
Period  of  oscillation 
Absolute  temperature 

Period  of  time  during  which  intensity  of  wind  stress 
increases  linearly  across  the  stress  band 

Velocity  in  x direction 

Mean  wind  speed 

Depth-dependent  perturbations  in  horizontal  velocity  in 
x direction 

x-component  of  depth-averaged,  horizontal  velocity 
Layer-averaged  u-component  of  horizontal  velocity 
Characteristic  velocity  in  Lake  Erie 
Velocity  in  y direction 


A3 


t 


V 


n 


V 

s 

w 


w 

a 


x,y,z 

[z] 

al,a2,a3 

r 

c 

c 

^R 

n(x,y,t) 


X' 


[A] 

ii 

v 
£ 
P 
P 
P ’ 


Depth-dependent  perturbations  in  horizontal  velocity  in 
y direction 

y-component  of  depth-averaged,  horizontal  velocity 
Velocity  normal  to  shoreline 
Constant  translation  speed  of  wind  field 
Velocity  in  z direction 
Reference  vertical  velocity 

Wind  velocity  at  a certain  height  above  the  water  surface 
Cartesian  coordinates 

Modal  matrix,  where  columns  are  eigenvectors  {Z.}  of  [H] 
Coefficients  of  Taylor's  expansion 

u :+  iv  , a complex  variable  representing  horizontal  velocity 
Surface  displacement  as  a function  of  x and  y 
Vertical  component  of  vorticity 
Reference  surface  displacement 

Surface  displacement  as  a function  of  time  and  space 

2 

An  eigenvalue  expressed  as  w /g;  ttH/D  a constant  which 
is  proportional  to  the  ratio  of  water  depth  to  frictional 
depth 

Function  of  both  temperature  and  salinity 

Diagonal  matrix  of  eigenvalues 

Absolute  viscosity 

Coefficient  of  interfacial  friction 

Wave  amplitude 

Density 

Constant  reference  density 

Departure  from  p depending  on  salinity 

AU 


i 


' 


t 


p 

p 


a 

k 


o 

a 


T 


T 


wx 


T 

,T 


Ey 

wR 

wy 

T 

y 


T 


XZ 


4> 

[*] 

i|*(x,y) 


<J> 

b) 


Density  of  the  air;  atmospheric  density 
Depth-averaged  density  for  the  k^'  layer 
Transformed  vertical  coordinate 
Bottom  friction  coefficient 

x-  and  y-  components,  respectively,  of  bottom  friction 
Reference  wind  shear  stress 

x-  and  y-  components,  respectively,  of  wind  stress  at  the 
water  surface 

x and  y components  of  shear  stress 

x-  and  z-  components  of  shear  stress  interfacial 

Functional  integral 

Volume  transport  stream  function 

Vector  column  matrix  for  the  stream  function 

Mass  flux  stream  function 

Angle  between  wind  direction  and  y-axis 

Heat  dissipation 

2 ir/T 

Points  at  which  variables  must  be  interpreted  since  they 
are  not  in  the  staggered  grid 
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In  accordance  with  1®  70-2-3,  paragraph  6c(l)(b), 
dated  15  February  1973.  a facsimile  catalog  card 
in  Library  of  Congress  format  is  reproduced  below. 


Raney,  Donald  C 

Lake  Erie  International  Jetport  model  feasibility  in- 
vestigation; Report  17-4:  Numerical  model  feasibility 
study,  by  Donald  C.  Raney,  Donald  L.  Durham,  and 
H.  Lee  Butler.  Vicksburg,  U.  S.  Army  Engineer  Waterways 
Experiment  Station,  1977. 

1 v.  (various  pagings)  illus.  27  cm.  (U.  S. 
Waterways  Experiment  Station.  Technical  report  H-74-6, 
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